* Bug 15662 fixed: display in mode format(e) repaired
[scilab.git] / scilab / modules / ast / src / c / types / dtoa.c
1 #define IEEE_8087
2
3 /****************************************************************
4  *
5  * The author of this software is David M. Gay.
6  *
7  * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
8  *
9  * Permission to use, copy, modify, and distribute this software for any
10  * purpose without fee is hereby granted, provided that this entire notice
11  * is included in all copies of any software which is or includes a copy
12  * or modification of this software and in all copies of the supporting
13  * documentation for such software.
14  *
15  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
16  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
17  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
18  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
19  *
20  ***************************************************************/
21
22 /* Please send bug reports to David M. Gay (dmg at acm dot org,
23  * with " at " changed at "@" and " dot " changed to ".").      */
24
25 /* On a machine with IEEE extended-precision registers, it is
26  * necessary to specify double-precision (53-bit) rounding precision
27  * before invoking strtod or dtoa.  If the machine uses (the equivalent
28  * of) Intel 80x87 arithmetic, the call
29  *      _control87(PC_53, MCW_PC);
30  * does this with many compilers.  Whether this or another call is
31  * appropriate depends on the compiler; for this to work, it may be
32  * necessary to #include "float.h" or another system-dependent header
33  * file.
34  */
35
36 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
37  * (Note that IEEE arithmetic is disabled by gcc's -ffast-math flag.)
38  *
39  * This strtod returns a nearest machine number to the input decimal
40  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
41  * broken by the IEEE round-even rule.  Otherwise ties are broken by
42  * biased rounding (add half and chop).
43  *
44  * Inspired loosely by William D. Clinger's paper "How to Read Floating
45  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
46  *
47  * Modifications:
48  *
49  *      1. We only require IEEE, IBM, or VAX double-precision
50  *              arithmetic (not IEEE double-extended).
51  *      2. We get by with floating-point arithmetic in a case that
52  *              Clinger missed -- when we're computing d * 10^n
53  *              for a small integer d and the integer n is not too
54  *              much larger than 22 (the maximum integer k for which
55  *              we can represent 10^k exactly), we may be able to
56  *              compute (d*10^k) * 10^(e-k) with just one roundoff.
57  *      3. Rather than a bit-at-a-time adjustment of the binary
58  *              result in the hard case, we use floating-point
59  *              arithmetic to determine the adjustment to within
60  *              one bit; only in really hard cases do we need to
61  *              compute a second residual.
62  *      4. Because of 3., we don't need a large table of powers of 10
63  *              for ten-to-e (just some small tables, e.g. of 10^k
64  *              for 0 <= k <= 22).
65  */
66
67 /*
68  * #define IEEE_8087 for IEEE-arithmetic machines where the least
69  *      significant byte has the lowest address.
70  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
71  *      significant byte has the lowest address.
72  * #define Long int on machines with 32-bit ints and 64-bit longs.
73  * #define IBM for IBM mainframe-style floating-point arithmetic.
74  * #define VAX for VAX-style floating-point arithmetic (D_floating).
75  * #define No_leftright to omit left-right logic in fast floating-point
76  *      computation of dtoa.  This will cause dtoa modes 4 and 5 to be
77  *      treated the same as modes 2 and 3 for some inputs.
78  * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
79  *      and strtod and dtoa should round accordingly.  Unless Trust_FLT_ROUNDS
80  *      is also #defined, fegetround() will be queried for the rounding mode.
81  *      Note that both FLT_ROUNDS and fegetround() are specified by the C99
82  *      standard (and are specified to be consistent, with fesetround()
83  *      affecting the value of FLT_ROUNDS), but that some (Linux) systems
84  *      do not work correctly in this regard, so using fegetround() is more
85  *      portable than using FLT_ROUNDS directly.
86  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
87  *      and Honor_FLT_ROUNDS is not #defined.
88  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
89  *      that use extended-precision instructions to compute rounded
90  *      products and quotients) with IBM.
91  * #define ROUND_BIASED for IEEE-format with biased rounding and arithmetic
92  *      that rounds toward +Infinity.
93  * #define ROUND_BIASED_without_Round_Up for IEEE-format with biased
94  *      rounding when the underlying floating-point arithmetic uses
95  *      unbiased rounding.  This prevent using ordinary floating-point
96  *      arithmetic when the result could be computed with one rounding error.
97  * #define Inaccurate_Divide for IEEE-format with correctly rounded
98  *      products but inaccurate quotients, e.g., for Intel i860.
99  * #define NO_LONG_LONG on machines that do not have a "long long"
100  *      integer type (of >= 64 bits).  On such machines, you can
101  *      #define Just_16 to store 16 bits per 32-bit Long when doing
102  *      high-precision integer arithmetic.  Whether this speeds things
103  *      up or slows things down depends on the machine and the number
104  *      being converted.  If long long is available and the name is
105  *      something other than "long long", #define Llong to be the name,
106  *      and if "unsigned Llong" does not work as an unsigned version of
107  *      Llong, #define #ULLong to be the corresponding unsigned type.
108  * #define Bad_float_h if your system lacks a float.h or if it does not
109  *      define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
110  *      FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
111  * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
112  *      if memory is available and otherwise does something you deem
113  *      appropriate.  If MALLOC is undefined, malloc will be invoked
114  *      directly -- and assumed always to succeed.  Similarly, if you
115  *      want something other than the system's free() to be called to
116  *      recycle memory acquired from MALLOC, #define FREE to be the
117  *      name of the alternate routine.  (FREE or free is only called in
118  *      pathological cases, e.g., in a dtoa call after a dtoa return in
119  *      mode 3 with thousands of digits requested.)
120  * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
121  *      memory allocations from a private pool of memory when possible.
122  *      When used, the private pool is PRIVATE_MEM bytes long:  2304 bytes,
123  *      unless #defined to be a different length.  This default length
124  *      suffices to get rid of MALLOC calls except for unusual cases,
125  *      such as decimal-to-binary conversion of a very long string of
126  *      digits.  The longest string dtoa can return is about 751 bytes
127  *      long.  For conversions by strtod of strings of 800 digits and
128  *      all dtoa conversions in single-threaded executions with 8-byte
129  *      pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
130  *      pointers, PRIVATE_MEM >= 7112 appears adequate.
131  * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
132  *      #defined automatically on IEEE systems.  On such systems,
133  *      when INFNAN_CHECK is #defined, strtod checks
134  *      for Infinity and NaN (case insensitively).  On some systems
135  *      (e.g., some HP systems), it may be necessary to #define NAN_WORD0
136  *      appropriately -- to the most significant word of a quiet NaN.
137  *      (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
138  *      When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
139  *      strtod also accepts (case insensitively) strings of the form
140  *      NaN(x), where x is a string of hexadecimal digits and spaces;
141  *      if there is only one string of hexadecimal digits, it is taken
142  *      for the 52 fraction bits of the resulting NaN; if there are two
143  *      or more strings of hex digits, the first is for the high 20 bits,
144  *      the second and subsequent for the low 32 bits, with intervening
145  *      white space ignored; but if this results in none of the 52
146  *      fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
147  *      and NAN_WORD1 are used instead.
148  * #define MULTIPLE_THREADS if the system offers preemptively scheduled
149  *      multiple threads.  In this case, you must provide (or suitably
150  *      #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
151  *      by FREE_DTOA_LOCK(n) for n = 0 or 1.  (The second lock, accessed
152  *      in pow5mult, ensures lazy evaluation of only one copy of high
153  *      powers of 5; omitting this lock would introduce a small
154  *      probability of wasting memory, but would otherwise be harmless.)
155  *      You must also invoke freedtoa(s) to free the value s returned by
156  *      dtoa.  You may do so whether or not MULTIPLE_THREADS is #defined.
157
158  *      When MULTIPLE_THREADS is #defined, this source file provides
159  *              void set_max_dtoa_threads(unsigned int n);
160  *      and expects
161  *              unsigned int dtoa_get_threadno(void);
162  *      to be available (possibly provided by
163  *              #define dtoa_get_threadno omp_get_thread_num
164  *      if OpenMP is in use or by
165  *              #define dtoa_get_threadno pthread_self
166  *      if Pthreads is in use), to return the current thread number.
167  *      If set_max_dtoa_threads(n) was called and the current thread
168  *      number is k with k < n, then calls on ACQUIRE_DTOA_LOCK(...) and
169  *      FREE_DTOA_LOCK(...) are avoided; instead each thread with thread
170  *      number < n has a separate copy of relevant data structures.
171  *      After set_max_dtoa_threads(n), a call set_max_dtoa_threads(m)
172  *      with m <= n has has no effect, but a call with m > n is honored.
173  *      Such a call invokes REALLOC (assumed to be "realloc" if REALLOC
174  *      is not #defined) to extend the size of the relevant array.
175
176  * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
177  *      avoids underflows on inputs whose result does not underflow.
178  *      If you #define NO_IEEE_Scale on a machine that uses IEEE-format
179  *      floating-point numbers and flushes underflows to zero rather
180  *      than implementing gradual underflow, then you must also #define
181  *      Sudden_Underflow.
182  * #define USE_LOCALE to use the current locale's decimal_point value.
183  * #define SET_INEXACT if IEEE arithmetic is being used and extra
184  *      computation should be done to set the inexact flag when the
185  *      result is inexact and avoid setting inexact when the result
186  *      is exact.  In this case, dtoa.c must be compiled in
187  *      an environment, perhaps provided by #include "dtoa.c" in a
188  *      suitable wrapper, that defines two functions,
189  *              int get_inexact(void);
190  *              void clear_inexact(void);
191  *      such that get_inexact() returns a nonzero value if the
192  *      inexact bit is already set, and clear_inexact() sets the
193  *      inexact bit to 0.  When SET_INEXACT is #defined, strtod
194  *      also does extra computations to set the underflow and overflow
195  *      flags when appropriate (i.e., when the result is tiny and
196  *      inexact or when it is a numeric value rounded to +-infinity).
197  * #define NO_ERRNO if strtod should not assign errno = ERANGE when
198  *      the result overflows to +-Infinity or underflows to 0.
199  *      When errno should be assigned, under seemingly rare conditions
200  *      it may be necessary to define Set_errno(x) suitably, e.g., in
201  *      a local errno.h, such as
202  *              #include <errno.h>
203  *              #define Set_errno(x) _set_errno(x)
204  * #define NO_HEX_FP to omit recognition of hexadecimal floating-point
205  *      values by strtod.
206  * #define NO_STRTOD_BIGCOMP (on IEEE-arithmetic systems only for now)
207  *      to disable logic for "fast" testing of very long input strings
208  *      to strtod.  This testing proceeds by initially truncating the
209  *      input string, then if necessary comparing the whole string with
210  *      a decimal expansion to decide close cases. This logic is only
211  *      used for input more than STRTOD_DIGLIM digits long (default 40).
212  */
213
214 #ifndef Long
215 #define Long int
216 #endif
217 #ifndef ULong
218 typedef unsigned Long ULong;
219 #endif
220
221 #ifdef DEBUG
222 #include <assert.h>
223 #include "stdio.h"
224 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
225 #define Debug(x) x
226 int dtoa_stats[7]; /* strtod_{64,96,bigcomp},dtoa_{exact,64,96,bigcomp} */
227 #else
228 #define assert(x) /*nothing*/
229 #define Debug(x) /*nothing*/
230 #endif
231
232 #include "stdlib.h"
233 #include "string.h"
234
235 #ifdef USE_LOCALE
236 #include "locale.h"
237 #endif
238
239 #ifdef Honor_FLT_ROUNDS
240 #ifndef Trust_FLT_ROUNDS
241 #include <fenv.h>
242 #endif
243 #endif
244
245 #ifdef __cplusplus
246 extern "C" {
247 #endif
248 #ifdef MALLOC
249 extern void *MALLOC(size_t);
250 #else
251 #define MALLOC malloc
252 #endif
253
254 #ifdef REALLOC
255 extern void *REALLOC(void*,size_t);
256 #else
257 #define REALLOC realloc
258 #endif
259
260 #ifndef FREE
261 #define FREE free
262 #endif
263
264 #ifdef __cplusplus
265         }
266 #endif
267
268 #ifndef Omit_Private_Memory
269 #ifndef PRIVATE_MEM
270 #define PRIVATE_MEM 2304
271 #endif
272 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
273 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
274 #endif
275
276 #undef IEEE_Arith
277 #undef Avoid_Underflow
278 #ifdef IEEE_MC68k
279 #define IEEE_Arith
280 #endif
281 #ifdef IEEE_8087
282 #define IEEE_Arith
283 #endif
284
285 #ifdef IEEE_Arith
286 #ifndef NO_INFNAN_CHECK
287 #undef INFNAN_CHECK
288 #define INFNAN_CHECK
289 #endif
290 #else
291 #undef INFNAN_CHECK
292 #define NO_STRTOD_BIGCOMP
293 #endif
294
295 #include "errno.h"
296
297 #ifdef NO_ERRNO /*{*/
298 #undef Set_errno
299 #define Set_errno(x)
300 #else
301 #ifndef Set_errno
302 #define Set_errno(x) errno = x
303 #endif
304 #endif /*}*/
305
306 #ifdef Bad_float_h
307
308 #ifdef IEEE_Arith
309 #define DBL_DIG 15
310 #define DBL_MAX_10_EXP 308
311 #define DBL_MAX_EXP 1024
312 #define FLT_RADIX 2
313 #endif /*IEEE_Arith*/
314
315 #ifdef IBM
316 #define DBL_DIG 16
317 #define DBL_MAX_10_EXP 75
318 #define DBL_MAX_EXP 63
319 #define FLT_RADIX 16
320 #define DBL_MAX 7.2370055773322621e+75
321 #endif
322
323 #ifdef VAX
324 #define DBL_DIG 16
325 #define DBL_MAX_10_EXP 38
326 #define DBL_MAX_EXP 127
327 #define FLT_RADIX 2
328 #define DBL_MAX 1.7014118346046923e+38
329 #endif
330
331 #ifndef LONG_MAX
332 #define LONG_MAX 2147483647
333 #endif
334
335 #else /* ifndef Bad_float_h */
336 #include "float.h"
337 #endif /* Bad_float_h */
338
339 #ifndef __MATH_H__
340 #include "math.h"
341 #endif
342
343 #ifdef __cplusplus
344 extern "C" {
345 #endif
346
347 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
348 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
349 #endif
350
351 #undef USE_BF96
352
353 #ifdef NO_LONG_LONG /*{{*/
354 #undef ULLong
355 #ifdef Just_16
356 #undef Pack_32
357 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
358  * This makes some inner loops simpler and sometimes saves work
359  * during multiplications, but it often seems to make things slightly
360  * slower.  Hence the default is now to store 32 bits per Long.
361  */
362 #endif
363 #else   /*}{ long long available */
364 #ifndef Llong
365 #define Llong long long
366 #endif
367 #ifndef ULLong
368 #define ULLong unsigned Llong
369 #endif
370 #ifndef NO_BF96 /*{*/
371 #define USE_BF96
372
373 #ifdef SET_INEXACT
374 #define dtoa_divmax 27
375 #else
376 int dtoa_divmax = 2;    /* Permit experimenting: on some systems, 64-bit integer */
377                         /* division is slow enough that we may sometimes want to */
378                         /* avoid using it.   We assume (but do not check) that   */
379                         /* dtoa_divmax <= 27.*/
380 #endif
381
382 typedef struct BF96 {           /* Normalized 96-bit software floating point numbers */
383         unsigned int b0,b1,b2;  /* b0 = most significant, binary point just to its left */
384         int e;                  /* number represented = b * 2^e, with .5 <= b < 1 */
385         } BF96;
386
387  static BF96 pten[667] = {
388         { 0xeef453d6, 0x923bd65a, 0x113faa29, -1136 },
389         { 0x9558b466, 0x1b6565f8, 0x4ac7ca59, -1132 },
390         { 0xbaaee17f, 0xa23ebf76, 0x5d79bcf0, -1129 },
391         { 0xe95a99df, 0x8ace6f53, 0xf4d82c2c, -1126 },
392         { 0x91d8a02b, 0xb6c10594, 0x79071b9b, -1122 },
393         { 0xb64ec836, 0xa47146f9, 0x9748e282, -1119 },
394         { 0xe3e27a44, 0x4d8d98b7, 0xfd1b1b23, -1116 },
395         { 0x8e6d8c6a, 0xb0787f72, 0xfe30f0f5, -1112 },
396         { 0xb208ef85, 0x5c969f4f, 0xbdbd2d33, -1109 },
397         { 0xde8b2b66, 0xb3bc4723, 0xad2c7880, -1106 },
398         { 0x8b16fb20, 0x3055ac76, 0x4c3bcb50, -1102 },
399         { 0xaddcb9e8, 0x3c6b1793, 0xdf4abe24, -1099 },
400         { 0xd953e862, 0x4b85dd78, 0xd71d6dad, -1096 },
401         { 0x87d4713d, 0x6f33aa6b, 0x8672648c, -1092 },
402         { 0xa9c98d8c, 0xcb009506, 0x680efdaf, -1089 },
403         { 0xd43bf0ef, 0xfdc0ba48, 0x0212bd1b, -1086 },
404         { 0x84a57695, 0xfe98746d, 0x014bb630, -1082 },
405         { 0xa5ced43b, 0x7e3e9188, 0x419ea3bd, -1079 },
406         { 0xcf42894a, 0x5dce35ea, 0x52064cac, -1076 },
407         { 0x818995ce, 0x7aa0e1b2, 0x7343efeb, -1072 },
408         { 0xa1ebfb42, 0x19491a1f, 0x1014ebe6, -1069 },
409         { 0xca66fa12, 0x9f9b60a6, 0xd41a26e0, -1066 },
410         { 0xfd00b897, 0x478238d0, 0x8920b098, -1063 },
411         { 0x9e20735e, 0x8cb16382, 0x55b46e5f, -1059 },
412         { 0xc5a89036, 0x2fddbc62, 0xeb2189f7, -1056 },
413         { 0xf712b443, 0xbbd52b7b, 0xa5e9ec75, -1053 },
414         { 0x9a6bb0aa, 0x55653b2d, 0x47b233c9, -1049 },
415         { 0xc1069cd4, 0xeabe89f8, 0x999ec0bb, -1046 },
416         { 0xf148440a, 0x256e2c76, 0xc00670ea, -1043 },
417         { 0x96cd2a86, 0x5764dbca, 0x38040692, -1039 },
418         { 0xbc807527, 0xed3e12bc, 0xc6050837, -1036 },
419         { 0xeba09271, 0xe88d976b, 0xf7864a44, -1033 },
420         { 0x93445b87, 0x31587ea3, 0x7ab3ee6a, -1029 },
421         { 0xb8157268, 0xfdae9e4c, 0x5960ea05, -1026 },
422         { 0xe61acf03, 0x3d1a45df, 0x6fb92487, -1023 },
423         { 0x8fd0c162, 0x06306bab, 0xa5d3b6d4, -1019 },
424         { 0xb3c4f1ba, 0x87bc8696, 0x8f48a489, -1016 },
425         { 0xe0b62e29, 0x29aba83c, 0x331acdab, -1013 },
426         { 0x8c71dcd9, 0xba0b4925, 0x9ff0c08b, -1009 },
427         { 0xaf8e5410, 0x288e1b6f, 0x07ecf0ae, -1006 },
428         { 0xdb71e914, 0x32b1a24a, 0xc9e82cd9, -1003 },
429         { 0x892731ac, 0x9faf056e, 0xbe311c08,  -999 },
430         { 0xab70fe17, 0xc79ac6ca, 0x6dbd630a,  -996 },
431         { 0xd64d3d9d, 0xb981787d, 0x092cbbcc,  -993 },
432         { 0x85f04682, 0x93f0eb4e, 0x25bbf560,  -989 },
433         { 0xa76c5823, 0x38ed2621, 0xaf2af2b8,  -986 },
434         { 0xd1476e2c, 0x07286faa, 0x1af5af66,  -983 },
435         { 0x82cca4db, 0x847945ca, 0x50d98d9f,  -979 },
436         { 0xa37fce12, 0x6597973c, 0xe50ff107,  -976 },
437         { 0xcc5fc196, 0xfefd7d0c, 0x1e53ed49,  -973 },
438         { 0xff77b1fc, 0xbebcdc4f, 0x25e8e89c,  -970 },
439         { 0x9faacf3d, 0xf73609b1, 0x77b19161,  -966 },
440         { 0xc795830d, 0x75038c1d, 0xd59df5b9,  -963 },
441         { 0xf97ae3d0, 0xd2446f25, 0x4b057328,  -960 },
442         { 0x9becce62, 0x836ac577, 0x4ee367f9,  -956 },
443         { 0xc2e801fb, 0x244576d5, 0x229c41f7,  -953 },
444         { 0xf3a20279, 0xed56d48a, 0x6b435275,  -950 },
445         { 0x9845418c, 0x345644d6, 0x830a1389,  -946 },
446         { 0xbe5691ef, 0x416bd60c, 0x23cc986b,  -943 },
447         { 0xedec366b, 0x11c6cb8f, 0x2cbfbe86,  -940 },
448         { 0x94b3a202, 0xeb1c3f39, 0x7bf7d714,  -936 },
449         { 0xb9e08a83, 0xa5e34f07, 0xdaf5ccd9,  -933 },
450         { 0xe858ad24, 0x8f5c22c9, 0xd1b3400f,  -930 },
451         { 0x91376c36, 0xd99995be, 0x23100809,  -926 },
452         { 0xb5854744, 0x8ffffb2d, 0xabd40a0c,  -923 },
453         { 0xe2e69915, 0xb3fff9f9, 0x16c90c8f,  -920 },
454         { 0x8dd01fad, 0x907ffc3b, 0xae3da7d9,  -916 },
455         { 0xb1442798, 0xf49ffb4a, 0x99cd11cf,  -913 },
456         { 0xdd95317f, 0x31c7fa1d, 0x40405643,  -910 },
457         { 0x8a7d3eef, 0x7f1cfc52, 0x482835ea,  -906 },
458         { 0xad1c8eab, 0x5ee43b66, 0xda324365,  -903 },
459         { 0xd863b256, 0x369d4a40, 0x90bed43e,  -900 },
460         { 0x873e4f75, 0xe2224e68, 0x5a7744a6,  -896 },
461         { 0xa90de353, 0x5aaae202, 0x711515d0,  -893 },
462         { 0xd3515c28, 0x31559a83, 0x0d5a5b44,  -890 },
463         { 0x8412d999, 0x1ed58091, 0xe858790a,  -886 },
464         { 0xa5178fff, 0x668ae0b6, 0x626e974d,  -883 },
465         { 0xce5d73ff, 0x402d98e3, 0xfb0a3d21,  -880 },
466         { 0x80fa687f, 0x881c7f8e, 0x7ce66634,  -876 },
467         { 0xa139029f, 0x6a239f72, 0x1c1fffc1,  -873 },
468         { 0xc9874347, 0x44ac874e, 0xa327ffb2,  -870 },
469         { 0xfbe91419, 0x15d7a922, 0x4bf1ff9f,  -867 },
470         { 0x9d71ac8f, 0xada6c9b5, 0x6f773fc3,  -863 },
471         { 0xc4ce17b3, 0x99107c22, 0xcb550fb4,  -860 },
472         { 0xf6019da0, 0x7f549b2b, 0x7e2a53a1,  -857 },
473         { 0x99c10284, 0x4f94e0fb, 0x2eda7444,  -853 },
474         { 0xc0314325, 0x637a1939, 0xfa911155,  -850 },
475         { 0xf03d93ee, 0xbc589f88, 0x793555ab,  -847 },
476         { 0x96267c75, 0x35b763b5, 0x4bc1558b,  -843 },
477         { 0xbbb01b92, 0x83253ca2, 0x9eb1aaed,  -840 },
478         { 0xea9c2277, 0x23ee8bcb, 0x465e15a9,  -837 },
479         { 0x92a1958a, 0x7675175f, 0x0bfacd89,  -833 },
480         { 0xb749faed, 0x14125d36, 0xcef980ec,  -830 },
481         { 0xe51c79a8, 0x5916f484, 0x82b7e127,  -827 },
482         { 0x8f31cc09, 0x37ae58d2, 0xd1b2ecb8,  -823 },
483         { 0xb2fe3f0b, 0x8599ef07, 0x861fa7e6,  -820 },
484         { 0xdfbdcece, 0x67006ac9, 0x67a791e0,  -817 },
485         { 0x8bd6a141, 0x006042bd, 0xe0c8bb2c,  -813 },
486         { 0xaecc4991, 0x4078536d, 0x58fae9f7,  -810 },
487         { 0xda7f5bf5, 0x90966848, 0xaf39a475,  -807 },
488         { 0x888f9979, 0x7a5e012d, 0x6d8406c9,  -803 },
489         { 0xaab37fd7, 0xd8f58178, 0xc8e5087b,  -800 },
490         { 0xd5605fcd, 0xcf32e1d6, 0xfb1e4a9a,  -797 },
491         { 0x855c3be0, 0xa17fcd26, 0x5cf2eea0,  -793 },
492         { 0xa6b34ad8, 0xc9dfc06f, 0xf42faa48,  -790 },
493         { 0xd0601d8e, 0xfc57b08b, 0xf13b94da,  -787 },
494         { 0x823c1279, 0x5db6ce57, 0x76c53d08,  -783 },
495         { 0xa2cb1717, 0xb52481ed, 0x54768c4b,  -780 },
496         { 0xcb7ddcdd, 0xa26da268, 0xa9942f5d,  -777 },
497         { 0xfe5d5415, 0x0b090b02, 0xd3f93b35,  -774 },
498         { 0x9efa548d, 0x26e5a6e1, 0xc47bc501,  -770 },
499         { 0xc6b8e9b0, 0x709f109a, 0x359ab641,  -767 },
500         { 0xf867241c, 0x8cc6d4c0, 0xc30163d2,  -764 },
501         { 0x9b407691, 0xd7fc44f8, 0x79e0de63,  -760 },
502         { 0xc2109436, 0x4dfb5636, 0x985915fc,  -757 },
503         { 0xf294b943, 0xe17a2bc4, 0x3e6f5b7b,  -754 },
504         { 0x979cf3ca, 0x6cec5b5a, 0xa705992c,  -750 },
505         { 0xbd8430bd, 0x08277231, 0x50c6ff78,  -747 },
506         { 0xece53cec, 0x4a314ebd, 0xa4f8bf56,  -744 },
507         { 0x940f4613, 0xae5ed136, 0x871b7795,  -740 },
508         { 0xb9131798, 0x99f68584, 0x28e2557b,  -737 },
509         { 0xe757dd7e, 0xc07426e5, 0x331aeada,  -734 },
510         { 0x9096ea6f, 0x3848984f, 0x3ff0d2c8,  -730 },
511         { 0xb4bca50b, 0x065abe63, 0x0fed077a,  -727 },
512         { 0xe1ebce4d, 0xc7f16dfb, 0xd3e84959,  -724 },
513         { 0x8d3360f0, 0x9cf6e4bd, 0x64712dd7,  -720 },
514         { 0xb080392c, 0xc4349dec, 0xbd8d794d,  -717 },
515         { 0xdca04777, 0xf541c567, 0xecf0d7a0,  -714 },
516         { 0x89e42caa, 0xf9491b60, 0xf41686c4,  -710 },
517         { 0xac5d37d5, 0xb79b6239, 0x311c2875,  -707 },
518         { 0xd77485cb, 0x25823ac7, 0x7d633293,  -704 },
519         { 0x86a8d39e, 0xf77164bc, 0xae5dff9c,  -700 },
520         { 0xa8530886, 0xb54dbdeb, 0xd9f57f83,  -697 },
521         { 0xd267caa8, 0x62a12d66, 0xd072df63,  -694 },
522         { 0x8380dea9, 0x3da4bc60, 0x4247cb9e,  -690 },
523         { 0xa4611653, 0x8d0deb78, 0x52d9be85,  -687 },
524         { 0xcd795be8, 0x70516656, 0x67902e27,  -684 },
525         { 0x806bd971, 0x4632dff6, 0x00ba1cd8,  -680 },
526         { 0xa086cfcd, 0x97bf97f3, 0x80e8a40e,  -677 },
527         { 0xc8a883c0, 0xfdaf7df0, 0x6122cd12,  -674 },
528         { 0xfad2a4b1, 0x3d1b5d6c, 0x796b8057,  -671 },
529         { 0x9cc3a6ee, 0xc6311a63, 0xcbe33036,  -667 },
530         { 0xc3f490aa, 0x77bd60fc, 0xbedbfc44,  -664 },
531         { 0xf4f1b4d5, 0x15acb93b, 0xee92fb55,  -661 },
532         { 0x99171105, 0x2d8bf3c5, 0x751bdd15,  -657 },
533         { 0xbf5cd546, 0x78eef0b6, 0xd262d45a,  -654 },
534         { 0xef340a98, 0x172aace4, 0x86fb8971,  -651 },
535         { 0x9580869f, 0x0e7aac0e, 0xd45d35e6,  -647 },
536         { 0xbae0a846, 0xd2195712, 0x89748360,  -644 },
537         { 0xe998d258, 0x869facd7, 0x2bd1a438,  -641 },
538         { 0x91ff8377, 0x5423cc06, 0x7b6306a3,  -637 },
539         { 0xb67f6455, 0x292cbf08, 0x1a3bc84c,  -634 },
540         { 0xe41f3d6a, 0x7377eeca, 0x20caba5f,  -631 },
541         { 0x8e938662, 0x882af53e, 0x547eb47b,  -627 },
542         { 0xb23867fb, 0x2a35b28d, 0xe99e619a,  -624 },
543         { 0xdec681f9, 0xf4c31f31, 0x6405fa00,  -621 },
544         { 0x8b3c113c, 0x38f9f37e, 0xde83bc40,  -617 },
545         { 0xae0b158b, 0x4738705e, 0x9624ab50,  -614 },
546         { 0xd98ddaee, 0x19068c76, 0x3badd624,  -611 },
547         { 0x87f8a8d4, 0xcfa417c9, 0xe54ca5d7,  -607 },
548         { 0xa9f6d30a, 0x038d1dbc, 0x5e9fcf4c,  -604 },
549         { 0xd47487cc, 0x8470652b, 0x7647c320,  -601 },
550         { 0x84c8d4df, 0xd2c63f3b, 0x29ecd9f4,  -597 },
551         { 0xa5fb0a17, 0xc777cf09, 0xf4681071,  -594 },
552         { 0xcf79cc9d, 0xb955c2cc, 0x7182148d,  -591 },
553         { 0x81ac1fe2, 0x93d599bf, 0xc6f14cd8,  -587 },
554         { 0xa21727db, 0x38cb002f, 0xb8ada00e,  -584 },
555         { 0xca9cf1d2, 0x06fdc03b, 0xa6d90811,  -581 },
556         { 0xfd442e46, 0x88bd304a, 0x908f4a16,  -578 },
557         { 0x9e4a9cec, 0x15763e2e, 0x9a598e4e,  -574 },
558         { 0xc5dd4427, 0x1ad3cdba, 0x40eff1e1,  -571 },
559         { 0xf7549530, 0xe188c128, 0xd12bee59,  -568 },
560         { 0x9a94dd3e, 0x8cf578b9, 0x82bb74f8,  -564 },
561         { 0xc13a148e, 0x3032d6e7, 0xe36a5236,  -561 },
562         { 0xf18899b1, 0xbc3f8ca1, 0xdc44e6c3,  -558 },
563         { 0x96f5600f, 0x15a7b7e5, 0x29ab103a,  -554 },
564         { 0xbcb2b812, 0xdb11a5de, 0x7415d448,  -551 },
565         { 0xebdf6617, 0x91d60f56, 0x111b495b,  -548 },
566         { 0x936b9fce, 0xbb25c995, 0xcab10dd9,  -544 },
567         { 0xb84687c2, 0x69ef3bfb, 0x3d5d514f,  -541 },
568         { 0xe65829b3, 0x046b0afa, 0x0cb4a5a3,  -538 },
569         { 0x8ff71a0f, 0xe2c2e6dc, 0x47f0e785,  -534 },
570         { 0xb3f4e093, 0xdb73a093, 0x59ed2167,  -531 },
571         { 0xe0f218b8, 0xd25088b8, 0x306869c1,  -528 },
572         { 0x8c974f73, 0x83725573, 0x1e414218,  -524 },
573         { 0xafbd2350, 0x644eeacf, 0xe5d1929e,  -521 },
574         { 0xdbac6c24, 0x7d62a583, 0xdf45f746,  -518 },
575         { 0x894bc396, 0xce5da772, 0x6b8bba8c,  -514 },
576         { 0xab9eb47c, 0x81f5114f, 0x066ea92f,  -511 },
577         { 0xd686619b, 0xa27255a2, 0xc80a537b,  -508 },
578         { 0x8613fd01, 0x45877585, 0xbd06742c,  -504 },
579         { 0xa798fc41, 0x96e952e7, 0x2c481138,  -501 },
580         { 0xd17f3b51, 0xfca3a7a0, 0xf75a1586,  -498 },
581         { 0x82ef8513, 0x3de648c4, 0x9a984d73,  -494 },
582         { 0xa3ab6658, 0x0d5fdaf5, 0xc13e60d0,  -491 },
583         { 0xcc963fee, 0x10b7d1b3, 0x318df905,  -488 },
584         { 0xffbbcfe9, 0x94e5c61f, 0xfdf17746,  -485 },
585         { 0x9fd561f1, 0xfd0f9bd3, 0xfeb6ea8b,  -481 },
586         { 0xc7caba6e, 0x7c5382c8, 0xfe64a52e,  -478 },
587         { 0xf9bd690a, 0x1b68637b, 0x3dfdce7a,  -475 },
588         { 0x9c1661a6, 0x51213e2d, 0x06bea10c,  -471 },
589         { 0xc31bfa0f, 0xe5698db8, 0x486e494f,  -468 },
590         { 0xf3e2f893, 0xdec3f126, 0x5a89dba3,  -465 },
591         { 0x986ddb5c, 0x6b3a76b7, 0xf8962946,  -461 },
592         { 0xbe895233, 0x86091465, 0xf6bbb397,  -458 },
593         { 0xee2ba6c0, 0x678b597f, 0x746aa07d,  -455 },
594         { 0x94db4838, 0x40b717ef, 0xa8c2a44e,  -451 },
595         { 0xba121a46, 0x50e4ddeb, 0x92f34d62,  -448 },
596         { 0xe896a0d7, 0xe51e1566, 0x77b020ba,  -445 },
597         { 0x915e2486, 0xef32cd60, 0x0ace1474,  -441 },
598         { 0xb5b5ada8, 0xaaff80b8, 0x0d819992,  -438 },
599         { 0xe3231912, 0xd5bf60e6, 0x10e1fff6,  -435 },
600         { 0x8df5efab, 0xc5979c8f, 0xca8d3ffa,  -431 },
601         { 0xb1736b96, 0xb6fd83b3, 0xbd308ff8,  -428 },
602         { 0xddd0467c, 0x64bce4a0, 0xac7cb3f6,  -425 },
603         { 0x8aa22c0d, 0xbef60ee4, 0x6bcdf07a,  -421 },
604         { 0xad4ab711, 0x2eb3929d, 0x86c16c98,  -418 },
605         { 0xd89d64d5, 0x7a607744, 0xe871c7bf,  -415 },
606         { 0x87625f05, 0x6c7c4a8b, 0x11471cd7,  -411 },
607         { 0xa93af6c6, 0xc79b5d2d, 0xd598e40d,  -408 },
608         { 0xd389b478, 0x79823479, 0x4aff1d10,  -405 },
609         { 0x843610cb, 0x4bf160cb, 0xcedf722a,  -401 },
610         { 0xa54394fe, 0x1eedb8fe, 0xc2974eb4,  -398 },
611         { 0xce947a3d, 0xa6a9273e, 0x733d2262,  -395 },
612         { 0x811ccc66, 0x8829b887, 0x0806357d,  -391 },
613         { 0xa163ff80, 0x2a3426a8, 0xca07c2dc,  -388 },
614         { 0xc9bcff60, 0x34c13052, 0xfc89b393,  -385 },
615         { 0xfc2c3f38, 0x41f17c67, 0xbbac2078,  -382 },
616         { 0x9d9ba783, 0x2936edc0, 0xd54b944b,  -378 },
617         { 0xc5029163, 0xf384a931, 0x0a9e795e,  -375 },
618         { 0xf64335bc, 0xf065d37d, 0x4d4617b5,  -372 },
619         { 0x99ea0196, 0x163fa42e, 0x504bced1,  -368 },
620         { 0xc06481fb, 0x9bcf8d39, 0xe45ec286,  -365 },
621         { 0xf07da27a, 0x82c37088, 0x5d767327,  -362 },
622         { 0x964e858c, 0x91ba2655, 0x3a6a07f8,  -358 },
623         { 0xbbe226ef, 0xb628afea, 0x890489f7,  -355 },
624         { 0xeadab0ab, 0xa3b2dbe5, 0x2b45ac74,  -352 },
625         { 0x92c8ae6b, 0x464fc96f, 0x3b0b8bc9,  -348 },
626         { 0xb77ada06, 0x17e3bbcb, 0x09ce6ebb,  -345 },
627         { 0xe5599087, 0x9ddcaabd, 0xcc420a6a,  -342 },
628         { 0x8f57fa54, 0xc2a9eab6, 0x9fa94682,  -338 },
629         { 0xb32df8e9, 0xf3546564, 0x47939822,  -335 },
630         { 0xdff97724, 0x70297ebd, 0x59787e2b,  -332 },
631         { 0x8bfbea76, 0xc619ef36, 0x57eb4edb,  -328 },
632         { 0xaefae514, 0x77a06b03, 0xede62292,  -325 },
633         { 0xdab99e59, 0x958885c4, 0xe95fab36,  -322 },
634         { 0x88b402f7, 0xfd75539b, 0x11dbcb02,  -318 },
635         { 0xaae103b5, 0xfcd2a881, 0xd652bdc2,  -315 },
636         { 0xd59944a3, 0x7c0752a2, 0x4be76d33,  -312 },
637         { 0x857fcae6, 0x2d8493a5, 0x6f70a440,  -308 },
638         { 0xa6dfbd9f, 0xb8e5b88e, 0xcb4ccd50,  -305 },
639         { 0xd097ad07, 0xa71f26b2, 0x7e2000a4,  -302 },
640         { 0x825ecc24, 0xc873782f, 0x8ed40066,  -298 },
641         { 0xa2f67f2d, 0xfa90563b, 0x72890080,  -295 },
642         { 0xcbb41ef9, 0x79346bca, 0x4f2b40a0,  -292 },
643         { 0xfea126b7, 0xd78186bc, 0xe2f610c8,  -289 },
644         { 0x9f24b832, 0xe6b0f436, 0x0dd9ca7d,  -285 },
645         { 0xc6ede63f, 0xa05d3143, 0x91503d1c,  -282 },
646         { 0xf8a95fcf, 0x88747d94, 0x75a44c63,  -279 },
647         { 0x9b69dbe1, 0xb548ce7c, 0xc986afbe,  -275 },
648         { 0xc24452da, 0x229b021b, 0xfbe85bad,  -272 },
649         { 0xf2d56790, 0xab41c2a2, 0xfae27299,  -269 },
650         { 0x97c560ba, 0x6b0919a5, 0xdccd879f,  -265 },
651         { 0xbdb6b8e9, 0x05cb600f, 0x5400e987,  -262 },
652         { 0xed246723, 0x473e3813, 0x290123e9,  -259 },
653         { 0x9436c076, 0x0c86e30b, 0xf9a0b672,  -255 },
654         { 0xb9447093, 0x8fa89bce, 0xf808e40e,  -252 },
655         { 0xe7958cb8, 0x7392c2c2, 0xb60b1d12,  -249 },
656         { 0x90bd77f3, 0x483bb9b9, 0xb1c6f22b,  -245 },
657         { 0xb4ecd5f0, 0x1a4aa828, 0x1e38aeb6,  -242 },
658         { 0xe2280b6c, 0x20dd5232, 0x25c6da63,  -239 },
659         { 0x8d590723, 0x948a535f, 0x579c487e,  -235 },
660         { 0xb0af48ec, 0x79ace837, 0x2d835a9d,  -232 },
661         { 0xdcdb1b27, 0x98182244, 0xf8e43145,  -229 },
662         { 0x8a08f0f8, 0xbf0f156b, 0x1b8e9ecb,  -225 },
663         { 0xac8b2d36, 0xeed2dac5, 0xe272467e,  -222 },
664         { 0xd7adf884, 0xaa879177, 0x5b0ed81d,  -219 },
665         { 0x86ccbb52, 0xea94baea, 0x98e94712,  -215 },
666         { 0xa87fea27, 0xa539e9a5, 0x3f2398d7,  -212 },
667         { 0xd29fe4b1, 0x8e88640e, 0x8eec7f0d,  -209 },
668         { 0x83a3eeee, 0xf9153e89, 0x1953cf68,  -205 },
669         { 0xa48ceaaa, 0xb75a8e2b, 0x5fa8c342,  -202 },
670         { 0xcdb02555, 0x653131b6, 0x3792f412,  -199 },
671         { 0x808e1755, 0x5f3ebf11, 0xe2bbd88b,  -195 },
672         { 0xa0b19d2a, 0xb70e6ed6, 0x5b6aceae,  -192 },
673         { 0xc8de0475, 0x64d20a8b, 0xf245825a,  -189 },
674         { 0xfb158592, 0xbe068d2e, 0xeed6e2f0,  -186 },
675         { 0x9ced737b, 0xb6c4183d, 0x55464dd6,  -182 },
676         { 0xc428d05a, 0xa4751e4c, 0xaa97e14c,  -179 },
677         { 0xf5330471, 0x4d9265df, 0xd53dd99f,  -176 },
678         { 0x993fe2c6, 0xd07b7fab, 0xe546a803,  -172 },
679         { 0xbf8fdb78, 0x849a5f96, 0xde985204,  -169 },
680         { 0xef73d256, 0xa5c0f77c, 0x963e6685,  -166 },
681         { 0x95a86376, 0x27989aad, 0xdde70013,  -162 },
682         { 0xbb127c53, 0xb17ec159, 0x5560c018,  -159 },
683         { 0xe9d71b68, 0x9dde71af, 0xaab8f01e,  -156 },
684         { 0x92267121, 0x62ab070d, 0xcab39613,  -152 },
685         { 0xb6b00d69, 0xbb55c8d1, 0x3d607b97,  -149 },
686         { 0xe45c10c4, 0x2a2b3b05, 0x8cb89a7d,  -146 },
687         { 0x8eb98a7a, 0x9a5b04e3, 0x77f3608e,  -142 },
688         { 0xb267ed19, 0x40f1c61c, 0x55f038b2,  -139 },
689         { 0xdf01e85f, 0x912e37a3, 0x6b6c46de,  -136 },
690         { 0x8b61313b, 0xbabce2c6, 0x2323ac4b,  -132 },
691         { 0xae397d8a, 0xa96c1b77, 0xabec975e,  -129 },
692         { 0xd9c7dced, 0x53c72255, 0x96e7bd35,  -126 },
693         { 0x881cea14, 0x545c7575, 0x7e50d641,  -122 },
694         { 0xaa242499, 0x697392d2, 0xdde50bd1,  -119 },
695         { 0xd4ad2dbf, 0xc3d07787, 0x955e4ec6,  -116 },
696         { 0x84ec3c97, 0xda624ab4, 0xbd5af13b,  -112 },
697         { 0xa6274bbd, 0xd0fadd61, 0xecb1ad8a,  -109 },
698         { 0xcfb11ead, 0x453994ba, 0x67de18ed,  -106 },
699         { 0x81ceb32c, 0x4b43fcf4, 0x80eacf94,  -102 },
700         { 0xa2425ff7, 0x5e14fc31, 0xa1258379,   -99 },
701         { 0xcad2f7f5, 0x359a3b3e, 0x096ee458,   -96 },
702         { 0xfd87b5f2, 0x8300ca0d, 0x8bca9d6e,   -93 },
703         { 0x9e74d1b7, 0x91e07e48, 0x775ea264,   -89 },
704         { 0xc6120625, 0x76589dda, 0x95364afe,   -86 },
705         { 0xf79687ae, 0xd3eec551, 0x3a83ddbd,   -83 },
706         { 0x9abe14cd, 0x44753b52, 0xc4926a96,   -79 },
707         { 0xc16d9a00, 0x95928a27, 0x75b7053c,   -76 },
708         { 0xf1c90080, 0xbaf72cb1, 0x5324c68b,   -73 },
709         { 0x971da050, 0x74da7bee, 0xd3f6fc16,   -69 },
710         { 0xbce50864, 0x92111aea, 0x88f4bb1c,   -66 },
711         { 0xec1e4a7d, 0xb69561a5, 0x2b31e9e3,   -63 },
712         { 0x9392ee8e, 0x921d5d07, 0x3aff322e,   -59 },
713         { 0xb877aa32, 0x36a4b449, 0x09befeb9,   -56 },
714         { 0xe69594be, 0xc44de15b, 0x4c2ebe68,   -53 },
715         { 0x901d7cf7, 0x3ab0acd9, 0x0f9d3701,   -49 },
716         { 0xb424dc35, 0x095cd80f, 0x538484c1,   -46 },
717         { 0xe12e1342, 0x4bb40e13, 0x2865a5f2,   -43 },
718         { 0x8cbccc09, 0x6f5088cb, 0xf93f87b7,   -39 },
719         { 0xafebff0b, 0xcb24aafe, 0xf78f69a5,   -36 },
720         { 0xdbe6fece, 0xbdedd5be, 0xb573440e,   -33 },
721         { 0x89705f41, 0x36b4a597, 0x31680a88,   -29 },
722         { 0xabcc7711, 0x8461cefc, 0xfdc20d2b,   -26 },
723         { 0xd6bf94d5, 0xe57a42bc, 0x3d329076,   -23 },
724         { 0x8637bd05, 0xaf6c69b5, 0xa63f9a49,   -19 },
725         { 0xa7c5ac47, 0x1b478423, 0x0fcf80dc,   -16 },
726         { 0xd1b71758, 0xe219652b, 0xd3c36113,   -13 },
727         { 0x83126e97, 0x8d4fdf3b, 0x645a1cac,    -9 },
728         { 0xa3d70a3d, 0x70a3d70a, 0x3d70a3d7,    -6 },
729         { 0xcccccccc, 0xcccccccc, 0xcccccccc,    -3 },
730         { 0x80000000, 0x00000000, 0x00000000,    1 },
731         { 0xa0000000, 0x00000000, 0x00000000,    4 },
732         { 0xc8000000, 0x00000000, 0x00000000,    7 },
733         { 0xfa000000, 0x00000000, 0x00000000,   10 },
734         { 0x9c400000, 0x00000000, 0x00000000,   14 },
735         { 0xc3500000, 0x00000000, 0x00000000,   17 },
736         { 0xf4240000, 0x00000000, 0x00000000,   20 },
737         { 0x98968000, 0x00000000, 0x00000000,   24 },
738         { 0xbebc2000, 0x00000000, 0x00000000,   27 },
739         { 0xee6b2800, 0x00000000, 0x00000000,   30 },
740         { 0x9502f900, 0x00000000, 0x00000000,   34 },
741         { 0xba43b740, 0x00000000, 0x00000000,   37 },
742         { 0xe8d4a510, 0x00000000, 0x00000000,   40 },
743         { 0x9184e72a, 0x00000000, 0x00000000,   44 },
744         { 0xb5e620f4, 0x80000000, 0x00000000,   47 },
745         { 0xe35fa931, 0xa0000000, 0x00000000,   50 },
746         { 0x8e1bc9bf, 0x04000000, 0x00000000,   54 },
747         { 0xb1a2bc2e, 0xc5000000, 0x00000000,   57 },
748         { 0xde0b6b3a, 0x76400000, 0x00000000,   60 },
749         { 0x8ac72304, 0x89e80000, 0x00000000,   64 },
750         { 0xad78ebc5, 0xac620000, 0x00000000,   67 },
751         { 0xd8d726b7, 0x177a8000, 0x00000000,   70 },
752         { 0x87867832, 0x6eac9000, 0x00000000,   74 },
753         { 0xa968163f, 0x0a57b400, 0x00000000,   77 },
754         { 0xd3c21bce, 0xcceda100, 0x00000000,   80 },
755         { 0x84595161, 0x401484a0, 0x00000000,   84 },
756         { 0xa56fa5b9, 0x9019a5c8, 0x00000000,   87 },
757         { 0xcecb8f27, 0xf4200f3a, 0x00000000,   90 },
758         { 0x813f3978, 0xf8940984, 0x40000000,   94 },
759         { 0xa18f07d7, 0x36b90be5, 0x50000000,   97 },
760         { 0xc9f2c9cd, 0x04674ede, 0xa4000000,  100 },
761         { 0xfc6f7c40, 0x45812296, 0x4d000000,  103 },
762         { 0x9dc5ada8, 0x2b70b59d, 0xf0200000,  107 },
763         { 0xc5371912, 0x364ce305, 0x6c280000,  110 },
764         { 0xf684df56, 0xc3e01bc6, 0xc7320000,  113 },
765         { 0x9a130b96, 0x3a6c115c, 0x3c7f4000,  117 },
766         { 0xc097ce7b, 0xc90715b3, 0x4b9f1000,  120 },
767         { 0xf0bdc21a, 0xbb48db20, 0x1e86d400,  123 },
768         { 0x96769950, 0xb50d88f4, 0x13144480,  127 },
769         { 0xbc143fa4, 0xe250eb31, 0x17d955a0,  130 },
770         { 0xeb194f8e, 0x1ae525fd, 0x5dcfab08,  133 },
771         { 0x92efd1b8, 0xd0cf37be, 0x5aa1cae5,  137 },
772         { 0xb7abc627, 0x050305ad, 0xf14a3d9e,  140 },
773         { 0xe596b7b0, 0xc643c719, 0x6d9ccd05,  143 },
774         { 0x8f7e32ce, 0x7bea5c6f, 0xe4820023,  147 },
775         { 0xb35dbf82, 0x1ae4f38b, 0xdda2802c,  150 },
776         { 0xe0352f62, 0xa19e306e, 0xd50b2037,  153 },
777         { 0x8c213d9d, 0xa502de45, 0x4526f422,  157 },
778         { 0xaf298d05, 0x0e4395d6, 0x9670b12b,  160 },
779         { 0xdaf3f046, 0x51d47b4c, 0x3c0cdd76,  163 },
780         { 0x88d8762b, 0xf324cd0f, 0xa5880a69,  167 },
781         { 0xab0e93b6, 0xefee0053, 0x8eea0d04,  170 },
782         { 0xd5d238a4, 0xabe98068, 0x72a49045,  173 },
783         { 0x85a36366, 0xeb71f041, 0x47a6da2b,  177 },
784         { 0xa70c3c40, 0xa64e6c51, 0x999090b6,  180 },
785         { 0xd0cf4b50, 0xcfe20765, 0xfff4b4e3,  183 },
786         { 0x82818f12, 0x81ed449f, 0xbff8f10e,  187 },
787         { 0xa321f2d7, 0x226895c7, 0xaff72d52,  190 },
788         { 0xcbea6f8c, 0xeb02bb39, 0x9bf4f8a6,  193 },
789         { 0xfee50b70, 0x25c36a08, 0x02f236d0,  196 },
790         { 0x9f4f2726, 0x179a2245, 0x01d76242,  200 },
791         { 0xc722f0ef, 0x9d80aad6, 0x424d3ad2,  203 },
792         { 0xf8ebad2b, 0x84e0d58b, 0xd2e08987,  206 },
793         { 0x9b934c3b, 0x330c8577, 0x63cc55f4,  210 },
794         { 0xc2781f49, 0xffcfa6d5, 0x3cbf6b71,  213 },
795         { 0xf316271c, 0x7fc3908a, 0x8bef464e,  216 },
796         { 0x97edd871, 0xcfda3a56, 0x97758bf0,  220 },
797         { 0xbde94e8e, 0x43d0c8ec, 0x3d52eeed,  223 },
798         { 0xed63a231, 0xd4c4fb27, 0x4ca7aaa8,  226 },
799         { 0x945e455f, 0x24fb1cf8, 0x8fe8caa9,  230 },
800         { 0xb975d6b6, 0xee39e436, 0xb3e2fd53,  233 },
801         { 0xe7d34c64, 0xa9c85d44, 0x60dbbca8,  236 },
802         { 0x90e40fbe, 0xea1d3a4a, 0xbc8955e9,  240 },
803         { 0xb51d13ae, 0xa4a488dd, 0x6babab63,  243 },
804         { 0xe264589a, 0x4dcdab14, 0xc696963c,  246 },
805         { 0x8d7eb760, 0x70a08aec, 0xfc1e1de5,  250 },
806         { 0xb0de6538, 0x8cc8ada8, 0x3b25a55f,  253 },
807         { 0xdd15fe86, 0xaffad912, 0x49ef0eb7,  256 },
808         { 0x8a2dbf14, 0x2dfcc7ab, 0x6e356932,  260 },
809         { 0xacb92ed9, 0x397bf996, 0x49c2c37f,  263 },
810         { 0xd7e77a8f, 0x87daf7fb, 0xdc33745e,  266 },
811         { 0x86f0ac99, 0xb4e8dafd, 0x69a028bb,  270 },
812         { 0xa8acd7c0, 0x222311bc, 0xc40832ea,  273 },
813         { 0xd2d80db0, 0x2aabd62b, 0xf50a3fa4,  276 },
814         { 0x83c7088e, 0x1aab65db, 0x792667c6,  280 },
815         { 0xa4b8cab1, 0xa1563f52, 0x577001b8,  283 },
816         { 0xcde6fd5e, 0x09abcf26, 0xed4c0226,  286 },
817         { 0x80b05e5a, 0xc60b6178, 0x544f8158,  290 },
818         { 0xa0dc75f1, 0x778e39d6, 0x696361ae,  293 },
819         { 0xc913936d, 0xd571c84c, 0x03bc3a19,  296 },
820         { 0xfb587849, 0x4ace3a5f, 0x04ab48a0,  299 },
821         { 0x9d174b2d, 0xcec0e47b, 0x62eb0d64,  303 },
822         { 0xc45d1df9, 0x42711d9a, 0x3ba5d0bd,  306 },
823         { 0xf5746577, 0x930d6500, 0xca8f44ec,  309 },
824         { 0x9968bf6a, 0xbbe85f20, 0x7e998b13,  313 },
825         { 0xbfc2ef45, 0x6ae276e8, 0x9e3fedd8,  316 },
826         { 0xefb3ab16, 0xc59b14a2, 0xc5cfe94e,  319 },
827         { 0x95d04aee, 0x3b80ece5, 0xbba1f1d1,  323 },
828         { 0xbb445da9, 0xca61281f, 0x2a8a6e45,  326 },
829         { 0xea157514, 0x3cf97226, 0xf52d09d7,  329 },
830         { 0x924d692c, 0xa61be758, 0x593c2626,  333 },
831         { 0xb6e0c377, 0xcfa2e12e, 0x6f8b2fb0,  336 },
832         { 0xe498f455, 0xc38b997a, 0x0b6dfb9c,  339 },
833         { 0x8edf98b5, 0x9a373fec, 0x4724bd41,  343 },
834         { 0xb2977ee3, 0x00c50fe7, 0x58edec91,  346 },
835         { 0xdf3d5e9b, 0xc0f653e1, 0x2f2967b6,  349 },
836         { 0x8b865b21, 0x5899f46c, 0xbd79e0d2,  353 },
837         { 0xae67f1e9, 0xaec07187, 0xecd85906,  356 },
838         { 0xda01ee64, 0x1a708de9, 0xe80e6f48,  359 },
839         { 0x884134fe, 0x908658b2, 0x3109058d,  363 },
840         { 0xaa51823e, 0x34a7eede, 0xbd4b46f0,  366 },
841         { 0xd4e5e2cd, 0xc1d1ea96, 0x6c9e18ac,  369 },
842         { 0x850fadc0, 0x9923329e, 0x03e2cf6b,  373 },
843         { 0xa6539930, 0xbf6bff45, 0x84db8346,  376 },
844         { 0xcfe87f7c, 0xef46ff16, 0xe6126418,  379 },
845         { 0x81f14fae, 0x158c5f6e, 0x4fcb7e8f,  383 },
846         { 0xa26da399, 0x9aef7749, 0xe3be5e33,  386 },
847         { 0xcb090c80, 0x01ab551c, 0x5cadf5bf,  389 },
848         { 0xfdcb4fa0, 0x02162a63, 0x73d9732f,  392 },
849         { 0x9e9f11c4, 0x014dda7e, 0x2867e7fd,  396 },
850         { 0xc646d635, 0x01a1511d, 0xb281e1fd,  399 },
851         { 0xf7d88bc2, 0x4209a565, 0x1f225a7c,  402 },
852         { 0x9ae75759, 0x6946075f, 0x3375788d,  406 },
853         { 0xc1a12d2f, 0xc3978937, 0x0052d6b1,  409 },
854         { 0xf209787b, 0xb47d6b84, 0xc0678c5d,  412 },
855         { 0x9745eb4d, 0x50ce6332, 0xf840b7ba,  416 },
856         { 0xbd176620, 0xa501fbff, 0xb650e5a9,  419 },
857         { 0xec5d3fa8, 0xce427aff, 0xa3e51f13,  422 },
858         { 0x93ba47c9, 0x80e98cdf, 0xc66f336c,  426 },
859         { 0xb8a8d9bb, 0xe123f017, 0xb80b0047,  429 },
860         { 0xe6d3102a, 0xd96cec1d, 0xa60dc059,  432 },
861         { 0x9043ea1a, 0xc7e41392, 0x87c89837,  436 },
862         { 0xb454e4a1, 0x79dd1877, 0x29babe45,  439 },
863         { 0xe16a1dc9, 0xd8545e94, 0xf4296dd6,  442 },
864         { 0x8ce2529e, 0x2734bb1d, 0x1899e4a6,  446 },
865         { 0xb01ae745, 0xb101e9e4, 0x5ec05dcf,  449 },
866         { 0xdc21a117, 0x1d42645d, 0x76707543,  452 },
867         { 0x899504ae, 0x72497eba, 0x6a06494a,  456 },
868         { 0xabfa45da, 0x0edbde69, 0x0487db9d,  459 },
869         { 0xd6f8d750, 0x9292d603, 0x45a9d284,  462 },
870         { 0x865b8692, 0x5b9bc5c2, 0x0b8a2392,  466 },
871         { 0xa7f26836, 0xf282b732, 0x8e6cac77,  469 },
872         { 0xd1ef0244, 0xaf2364ff, 0x3207d795,  472 },
873         { 0x8335616a, 0xed761f1f, 0x7f44e6bd,  476 },
874         { 0xa402b9c5, 0xa8d3a6e7, 0x5f16206c,  479 },
875         { 0xcd036837, 0x130890a1, 0x36dba887,  482 },
876         { 0x80222122, 0x6be55a64, 0xc2494954,  486 },
877         { 0xa02aa96b, 0x06deb0fd, 0xf2db9baa,  489 },
878         { 0xc83553c5, 0xc8965d3d, 0x6f928294,  492 },
879         { 0xfa42a8b7, 0x3abbf48c, 0xcb772339,  495 },
880         { 0x9c69a972, 0x84b578d7, 0xff2a7604,  499 },
881         { 0xc38413cf, 0x25e2d70d, 0xfef51385,  502 },
882         { 0xf46518c2, 0xef5b8cd1, 0x7eb25866,  505 },
883         { 0x98bf2f79, 0xd5993802, 0xef2f773f,  509 },
884         { 0xbeeefb58, 0x4aff8603, 0xaafb550f,  512 },
885         { 0xeeaaba2e, 0x5dbf6784, 0x95ba2a53,  515 },
886         { 0x952ab45c, 0xfa97a0b2, 0xdd945a74,  519 },
887         { 0xba756174, 0x393d88df, 0x94f97111,  522 },
888         { 0xe912b9d1, 0x478ceb17, 0x7a37cd56,  525 },
889         { 0x91abb422, 0xccb812ee, 0xac62e055,  529 },
890         { 0xb616a12b, 0x7fe617aa, 0x577b986b,  532 },
891         { 0xe39c4976, 0x5fdf9d94, 0xed5a7e85,  535 },
892         { 0x8e41ade9, 0xfbebc27d, 0x14588f13,  539 },
893         { 0xb1d21964, 0x7ae6b31c, 0x596eb2d8,  542 },
894         { 0xde469fbd, 0x99a05fe3, 0x6fca5f8e,  545 },
895         { 0x8aec23d6, 0x80043bee, 0x25de7bb9,  549 },
896         { 0xada72ccc, 0x20054ae9, 0xaf561aa7,  552 },
897         { 0xd910f7ff, 0x28069da4, 0x1b2ba151,  555 },
898         { 0x87aa9aff, 0x79042286, 0x90fb44d2,  559 },
899         { 0xa99541bf, 0x57452b28, 0x353a1607,  562 },
900         { 0xd3fa922f, 0x2d1675f2, 0x42889b89,  565 },
901         { 0x847c9b5d, 0x7c2e09b7, 0x69956135,  569 },
902         { 0xa59bc234, 0xdb398c25, 0x43fab983,  572 },
903         { 0xcf02b2c2, 0x1207ef2e, 0x94f967e4,  575 },
904         { 0x8161afb9, 0x4b44f57d, 0x1d1be0ee,  579 },
905         { 0xa1ba1ba7, 0x9e1632dc, 0x6462d92a,  582 },
906         { 0xca28a291, 0x859bbf93, 0x7d7b8f75,  585 },
907         { 0xfcb2cb35, 0xe702af78, 0x5cda7352,  588 },
908         { 0x9defbf01, 0xb061adab, 0x3a088813,  592 },
909         { 0xc56baec2, 0x1c7a1916, 0x088aaa18,  595 },
910         { 0xf6c69a72, 0xa3989f5b, 0x8aad549e,  598 },
911         { 0x9a3c2087, 0xa63f6399, 0x36ac54e2,  602 },
912         { 0xc0cb28a9, 0x8fcf3c7f, 0x84576a1b,  605 },
913         { 0xf0fdf2d3, 0xf3c30b9f, 0x656d44a2,  608 },
914         { 0x969eb7c4, 0x7859e743, 0x9f644ae5,  612 },
915         { 0xbc4665b5, 0x96706114, 0x873d5d9f,  615 },
916         { 0xeb57ff22, 0xfc0c7959, 0xa90cb506,  618 },
917         { 0x9316ff75, 0xdd87cbd8, 0x09a7f124,  622 },
918         { 0xb7dcbf53, 0x54e9bece, 0x0c11ed6d,  625 },
919         { 0xe5d3ef28, 0x2a242e81, 0x8f1668c8,  628 },
920         { 0x8fa47579, 0x1a569d10, 0xf96e017d,  632 },
921         { 0xb38d92d7, 0x60ec4455, 0x37c981dc,  635 },
922         { 0xe070f78d, 0x3927556a, 0x85bbe253,  638 },
923         { 0x8c469ab8, 0x43b89562, 0x93956d74,  642 },
924         { 0xaf584166, 0x54a6babb, 0x387ac8d1,  645 },
925         { 0xdb2e51bf, 0xe9d0696a, 0x06997b05,  648 },
926         { 0x88fcf317, 0xf22241e2, 0x441fece3,  652 },
927         { 0xab3c2fdd, 0xeeaad25a, 0xd527e81c,  655 },
928         { 0xd60b3bd5, 0x6a5586f1, 0x8a71e223,  658 },
929         { 0x85c70565, 0x62757456, 0xf6872d56,  662 },
930         { 0xa738c6be, 0xbb12d16c, 0xb428f8ac,  665 },
931         { 0xd106f86e, 0x69d785c7, 0xe13336d7,  668 },
932         { 0x82a45b45, 0x0226b39c, 0xecc00246,  672 },
933         { 0xa34d7216, 0x42b06084, 0x27f002d7,  675 },
934         { 0xcc20ce9b, 0xd35c78a5, 0x31ec038d,  678 },
935         { 0xff290242, 0xc83396ce, 0x7e670471,  681 },
936         { 0x9f79a169, 0xbd203e41, 0x0f0062c6,  685 },
937         { 0xc75809c4, 0x2c684dd1, 0x52c07b78,  688 },
938         { 0xf92e0c35, 0x37826145, 0xa7709a56,  691 },
939         { 0x9bbcc7a1, 0x42b17ccb, 0x88a66076,  695 },
940         { 0xc2abf989, 0x935ddbfe, 0x6acff893,  698 },
941         { 0xf356f7eb, 0xf83552fe, 0x0583f6b8,  701 },
942         { 0x98165af3, 0x7b2153de, 0xc3727a33,  705 },
943         { 0xbe1bf1b0, 0x59e9a8d6, 0x744f18c0,  708 },
944         { 0xeda2ee1c, 0x7064130c, 0x1162def0,  711 },
945         { 0x9485d4d1, 0xc63e8be7, 0x8addcb56,  715 },
946         { 0xb9a74a06, 0x37ce2ee1, 0x6d953e2b,  718 },
947         { 0xe8111c87, 0xc5c1ba99, 0xc8fa8db6,  721 },
948         { 0x910ab1d4, 0xdb9914a0, 0x1d9c9892,  725 },
949         { 0xb54d5e4a, 0x127f59c8, 0x2503beb6,  728 },
950         { 0xe2a0b5dc, 0x971f303a, 0x2e44ae64,  731 },
951         { 0x8da471a9, 0xde737e24, 0x5ceaecfe,  735 },
952         { 0xb10d8e14, 0x56105dad, 0x7425a83e,  738 },
953         { 0xdd50f199, 0x6b947518, 0xd12f124e,  741 },
954         { 0x8a5296ff, 0xe33cc92f, 0x82bd6b70,  745 },
955         { 0xace73cbf, 0xdc0bfb7b, 0x636cc64d,  748 },
956         { 0xd8210bef, 0xd30efa5a, 0x3c47f7e0,  751 },
957         { 0x8714a775, 0xe3e95c78, 0x65acfaec,  755 },
958         { 0xa8d9d153, 0x5ce3b396, 0x7f1839a7,  758 },
959         { 0xd31045a8, 0x341ca07c, 0x1ede4811,  761 },
960         { 0x83ea2b89, 0x2091e44d, 0x934aed0a,  765 },
961         { 0xa4e4b66b, 0x68b65d60, 0xf81da84d,  768 },
962         { 0xce1de406, 0x42e3f4b9, 0x36251260,  771 },
963         { 0x80d2ae83, 0xe9ce78f3, 0xc1d72b7c,  775 },
964         { 0xa1075a24, 0xe4421730, 0xb24cf65b,  778 },
965         { 0xc94930ae, 0x1d529cfc, 0xdee033f2,  781 },
966         { 0xfb9b7cd9, 0xa4a7443c, 0x169840ef,  784 },
967         { 0x9d412e08, 0x06e88aa5, 0x8e1f2895,  788 },
968         { 0xc491798a, 0x08a2ad4e, 0xf1a6f2ba,  791 },
969         { 0xf5b5d7ec, 0x8acb58a2, 0xae10af69,  794 },
970         { 0x9991a6f3, 0xd6bf1765, 0xacca6da1,  798 },
971         { 0xbff610b0, 0xcc6edd3f, 0x17fd090a,  801 },
972         { 0xeff394dc, 0xff8a948e, 0xddfc4b4c,  804 },
973         { 0x95f83d0a, 0x1fb69cd9, 0x4abdaf10,  808 },
974         { 0xbb764c4c, 0xa7a4440f, 0x9d6d1ad4,  811 },
975         { 0xea53df5f, 0xd18d5513, 0x84c86189,  814 },
976         { 0x92746b9b, 0xe2f8552c, 0x32fd3cf5,  818 },
977         { 0xb7118682, 0xdbb66a77, 0x3fbc8c33,  821 },
978         { 0xe4d5e823, 0x92a40515, 0x0fabaf3f,  824 },
979         { 0x8f05b116, 0x3ba6832d, 0x29cb4d87,  828 },
980         { 0xb2c71d5b, 0xca9023f8, 0x743e20e9,  831 },
981         { 0xdf78e4b2, 0xbd342cf6, 0x914da924,  834 },
982         { 0x8bab8eef, 0xb6409c1a, 0x1ad089b6,  838 },
983         { 0xae9672ab, 0xa3d0c320, 0xa184ac24,  841 },
984         { 0xda3c0f56, 0x8cc4f3e8, 0xc9e5d72d,  844 },
985         { 0x88658996, 0x17fb1871, 0x7e2fa67c,  848 },
986         { 0xaa7eebfb, 0x9df9de8d, 0xddbb901b,  851 },
987         { 0xd51ea6fa, 0x85785631, 0x552a7422,  854 },
988         { 0x8533285c, 0x936b35de, 0xd53a8895,  858 },
989         { 0xa67ff273, 0xb8460356, 0x8a892aba,  861 },
990         { 0xd01fef10, 0xa657842c, 0x2d2b7569,  864 },
991         { 0x8213f56a, 0x67f6b29b, 0x9c3b2962,  868 },
992         { 0xa298f2c5, 0x01f45f42, 0x8349f3ba,  871 },
993         { 0xcb3f2f76, 0x42717713, 0x241c70a9,  874 },
994         { 0xfe0efb53, 0xd30dd4d7, 0xed238cd3,  877 },
995         { 0x9ec95d14, 0x63e8a506, 0xf4363804,  881 },
996         { 0xc67bb459, 0x7ce2ce48, 0xb143c605,  884 },
997         { 0xf81aa16f, 0xdc1b81da, 0xdd94b786,  887 },
998         { 0x9b10a4e5, 0xe9913128, 0xca7cf2b4,  891 },
999         { 0xc1d4ce1f, 0x63f57d72, 0xfd1c2f61,  894 },
1000         { 0xf24a01a7, 0x3cf2dccf, 0xbc633b39,  897 },
1001         { 0x976e4108, 0x8617ca01, 0xd5be0503,  901 },
1002         { 0xbd49d14a, 0xa79dbc82, 0x4b2d8644,  904 },
1003         { 0xec9c459d, 0x51852ba2, 0xddf8e7d6,  907 },
1004         { 0x93e1ab82, 0x52f33b45, 0xcabb90e5,  911 },
1005         { 0xb8da1662, 0xe7b00a17, 0x3d6a751f,  914 },
1006         { 0xe7109bfb, 0xa19c0c9d, 0x0cc51267,  917 },
1007         { 0x906a617d, 0x450187e2, 0x27fb2b80,  921 },
1008         { 0xb484f9dc, 0x9641e9da, 0xb1f9f660,  924 },
1009         { 0xe1a63853, 0xbbd26451, 0x5e7873f8,  927 },
1010         { 0x8d07e334, 0x55637eb2, 0xdb0b487b,  931 },
1011         { 0xb049dc01, 0x6abc5e5f, 0x91ce1a9a,  934 },
1012         { 0xdc5c5301, 0xc56b75f7, 0x7641a140,  937 },
1013         { 0x89b9b3e1, 0x1b6329ba, 0xa9e904c8,  941 },
1014         { 0xac2820d9, 0x623bf429, 0x546345fa,  944 },
1015         { 0xd732290f, 0xbacaf133, 0xa97c1779,  947 },
1016         { 0x867f59a9, 0xd4bed6c0, 0x49ed8eab,  951 },
1017         { 0xa81f3014, 0x49ee8c70, 0x5c68f256,  954 },
1018         { 0xd226fc19, 0x5c6a2f8c, 0x73832eec,  957 },
1019         { 0x83585d8f, 0xd9c25db7, 0xc831fd53,  961 },
1020         { 0xa42e74f3, 0xd032f525, 0xba3e7ca8,  964 },
1021         { 0xcd3a1230, 0xc43fb26f, 0x28ce1bd2,  967 },
1022         { 0x80444b5e, 0x7aa7cf85, 0x7980d163,  971 },
1023         { 0xa0555e36, 0x1951c366, 0xd7e105bc,  974 },
1024         { 0xc86ab5c3, 0x9fa63440, 0x8dd9472b,  977 },
1025         { 0xfa856334, 0x878fc150, 0xb14f98f6,  980 },
1026         { 0x9c935e00, 0xd4b9d8d2, 0x6ed1bf9a,  984 },
1027         { 0xc3b83581, 0x09e84f07, 0x0a862f80,  987 },
1028         { 0xf4a642e1, 0x4c6262c8, 0xcd27bb61,  990 },
1029         { 0x98e7e9cc, 0xcfbd7dbd, 0x8038d51c,  994 },
1030         { 0xbf21e440, 0x03acdd2c, 0xe0470a63,  997 },
1031         { 0xeeea5d50, 0x04981478, 0x1858ccfc, 1000 },
1032         { 0x95527a52, 0x02df0ccb, 0x0f37801e, 1004 },
1033         { 0xbaa718e6, 0x8396cffd, 0xd3056025, 1007 },
1034         { 0xe950df20, 0x247c83fd, 0x47c6b82e, 1010 },
1035         { 0x91d28b74, 0x16cdd27e, 0x4cdc331d, 1014 },
1036         { 0xb6472e51, 0x1c81471d, 0xe0133fe4, 1017 },
1037         { 0xe3d8f9e5, 0x63a198e5, 0x58180fdd, 1020 },
1038         { 0x8e679c2f, 0x5e44ff8f, 0x570f09ea, 1024 },
1039         { 0xb201833b, 0x35d63f73, 0x2cd2cc65, 1027 },
1040         { 0xde81e40a, 0x034bcf4f, 0xf8077f7e, 1030 },
1041         { 0x8b112e86, 0x420f6191, 0xfb04afaf, 1034 },
1042         { 0xadd57a27, 0xd29339f6, 0x79c5db9a, 1037 },
1043         { 0xd94ad8b1, 0xc7380874, 0x18375281, 1040 },
1044         { 0x87cec76f, 0x1c830548, 0x8f229391, 1044 },
1045         { 0xa9c2794a, 0xe3a3c69a, 0xb2eb3875, 1047 },
1046         { 0xd433179d, 0x9c8cb841, 0x5fa60692, 1050 },
1047         { 0x849feec2, 0x81d7f328, 0xdbc7c41b, 1054 },
1048         { 0xa5c7ea73, 0x224deff3, 0x12b9b522, 1057 },
1049         { 0xcf39e50f, 0xeae16bef, 0xd768226b, 1060 },
1050         { 0x81842f29, 0xf2cce375, 0xe6a11583, 1064 },
1051         { 0xa1e53af4, 0x6f801c53, 0x60495ae3, 1067 },
1052         { 0xca5e89b1, 0x8b602368, 0x385bb19c, 1070 },
1053         { 0xfcf62c1d, 0xee382c42, 0x46729e03, 1073 },
1054         { 0x9e19db92, 0xb4e31ba9, 0x6c07a2c2, 1077 }
1055         };
1056  static short int Lhint[2098] = {
1057            /*18,*/19,    19,    19,    19,    20,    20,    20,    21,    21,
1058            21,    22,    22,    22,    23,    23,    23,    23,    24,    24,
1059            24,    25,    25,    25,    26,    26,    26,    26,    27,    27,
1060            27,    28,    28,    28,    29,    29,    29,    29,    30,    30,
1061            30,    31,    31,    31,    32,    32,    32,    32,    33,    33,
1062            33,    34,    34,    34,    35,    35,    35,    35,    36,    36,
1063            36,    37,    37,    37,    38,    38,    38,    38,    39,    39,
1064            39,    40,    40,    40,    41,    41,    41,    41,    42,    42,
1065            42,    43,    43,    43,    44,    44,    44,    44,    45,    45,
1066            45,    46,    46,    46,    47,    47,    47,    47,    48,    48,
1067            48,    49,    49,    49,    50,    50,    50,    51,    51,    51,
1068            51,    52,    52,    52,    53,    53,    53,    54,    54,    54,
1069            54,    55,    55,    55,    56,    56,    56,    57,    57,    57,
1070            57,    58,    58,    58,    59,    59,    59,    60,    60,    60,
1071            60,    61,    61,    61,    62,    62,    62,    63,    63,    63,
1072            63,    64,    64,    64,    65,    65,    65,    66,    66,    66,
1073            66,    67,    67,    67,    68,    68,    68,    69,    69,    69,
1074            69,    70,    70,    70,    71,    71,    71,    72,    72,    72,
1075            72,    73,    73,    73,    74,    74,    74,    75,    75,    75,
1076            75,    76,    76,    76,    77,    77,    77,    78,    78,    78,
1077            78,    79,    79,    79,    80,    80,    80,    81,    81,    81,
1078            82,    82,    82,    82,    83,    83,    83,    84,    84,    84,
1079            85,    85,    85,    85,    86,    86,    86,    87,    87,    87,
1080            88,    88,    88,    88,    89,    89,    89,    90,    90,    90,
1081            91,    91,    91,    91,    92,    92,    92,    93,    93,    93,
1082            94,    94,    94,    94,    95,    95,    95,    96,    96,    96,
1083            97,    97,    97,    97,    98,    98,    98,    99,    99,    99,
1084           100,   100,   100,   100,   101,   101,   101,   102,   102,   102,
1085           103,   103,   103,   103,   104,   104,   104,   105,   105,   105,
1086           106,   106,   106,   106,   107,   107,   107,   108,   108,   108,
1087           109,   109,   109,   110,   110,   110,   110,   111,   111,   111,
1088           112,   112,   112,   113,   113,   113,   113,   114,   114,   114,
1089           115,   115,   115,   116,   116,   116,   116,   117,   117,   117,
1090           118,   118,   118,   119,   119,   119,   119,   120,   120,   120,
1091           121,   121,   121,   122,   122,   122,   122,   123,   123,   123,
1092           124,   124,   124,   125,   125,   125,   125,   126,   126,   126,
1093           127,   127,   127,   128,   128,   128,   128,   129,   129,   129,
1094           130,   130,   130,   131,   131,   131,   131,   132,   132,   132,
1095           133,   133,   133,   134,   134,   134,   134,   135,   135,   135,
1096           136,   136,   136,   137,   137,   137,   137,   138,   138,   138,
1097           139,   139,   139,   140,   140,   140,   141,   141,   141,   141,
1098           142,   142,   142,   143,   143,   143,   144,   144,   144,   144,
1099           145,   145,   145,   146,   146,   146,   147,   147,   147,   147,
1100           148,   148,   148,   149,   149,   149,   150,   150,   150,   150,
1101           151,   151,   151,   152,   152,   152,   153,   153,   153,   153,
1102           154,   154,   154,   155,   155,   155,   156,   156,   156,   156,
1103           157,   157,   157,   158,   158,   158,   159,   159,   159,   159,
1104           160,   160,   160,   161,   161,   161,   162,   162,   162,   162,
1105           163,   163,   163,   164,   164,   164,   165,   165,   165,   165,
1106           166,   166,   166,   167,   167,   167,   168,   168,   168,   169,
1107           169,   169,   169,   170,   170,   170,   171,   171,   171,   172,
1108           172,   172,   172,   173,   173,   173,   174,   174,   174,   175,
1109           175,   175,   175,   176,   176,   176,   177,   177,   177,   178,
1110           178,   178,   178,   179,   179,   179,   180,   180,   180,   181,
1111           181,   181,   181,   182,   182,   182,   183,   183,   183,   184,
1112           184,   184,   184,   185,   185,   185,   186,   186,   186,   187,
1113           187,   187,   187,   188,   188,   188,   189,   189,   189,   190,
1114           190,   190,   190,   191,   191,   191,   192,   192,   192,   193,
1115           193,   193,   193,   194,   194,   194,   195,   195,   195,   196,
1116           196,   196,   197,   197,   197,   197,   198,   198,   198,   199,
1117           199,   199,   200,   200,   200,   200,   201,   201,   201,   202,
1118           202,   202,   203,   203,   203,   203,   204,   204,   204,   205,
1119           205,   205,   206,   206,   206,   206,   207,   207,   207,   208,
1120           208,   208,   209,   209,   209,   209,   210,   210,   210,   211,
1121           211,   211,   212,   212,   212,   212,   213,   213,   213,   214,
1122           214,   214,   215,   215,   215,   215,   216,   216,   216,   217,
1123           217,   217,   218,   218,   218,   218,   219,   219,   219,   220,
1124           220,   220,   221,   221,   221,   221,   222,   222,   222,   223,
1125           223,   223,   224,   224,   224,   224,   225,   225,   225,   226,
1126           226,   226,   227,   227,   227,   228,   228,   228,   228,   229,
1127           229,   229,   230,   230,   230,   231,   231,   231,   231,   232,
1128           232,   232,   233,   233,   233,   234,   234,   234,   234,   235,
1129           235,   235,   236,   236,   236,   237,   237,   237,   237,   238,
1130           238,   238,   239,   239,   239,   240,   240,   240,   240,   241,
1131           241,   241,   242,   242,   242,   243,   243,   243,   243,   244,
1132           244,   244,   245,   245,   245,   246,   246,   246,   246,   247,
1133           247,   247,   248,   248,   248,   249,   249,   249,   249,   250,
1134           250,   250,   251,   251,   251,   252,   252,   252,   252,   253,
1135           253,   253,   254,   254,   254,   255,   255,   255,   256,   256,
1136           256,   256,   257,   257,   257,   258,   258,   258,   259,   259,
1137           259,   259,   260,   260,   260,   261,   261,   261,   262,   262,
1138           262,   262,   263,   263,   263,   264,   264,   264,   265,   265,
1139           265,   265,   266,   266,   266,   267,   267,   267,   268,   268,
1140           268,   268,   269,   269,   269,   270,   270,   270,   271,   271,
1141           271,   271,   272,   272,   272,   273,   273,   273,   274,   274,
1142           274,   274,   275,   275,   275,   276,   276,   276,   277,   277,
1143           277,   277,   278,   278,   278,   279,   279,   279,   280,   280,
1144           280,   280,   281,   281,   281,   282,   282,   282,   283,   283,
1145           283,   283,   284,   284,   284,   285,   285,   285,   286,   286,
1146           286,   287,   287,   287,   287,   288,   288,   288,   289,   289,
1147           289,   290,   290,   290,   290,   291,   291,   291,   292,   292,
1148           292,   293,   293,   293,   293,   294,   294,   294,   295,   295,
1149           295,   296,   296,   296,   296,   297,   297,   297,   298,   298,
1150           298,   299,   299,   299,   299,   300,   300,   300,   301,   301,
1151           301,   302,   302,   302,   302,   303,   303,   303,   304,   304,
1152           304,   305,   305,   305,   305,   306,   306,   306,   307,   307,
1153           307,   308,   308,   308,   308,   309,   309,   309,   310,   310,
1154           310,   311,   311,   311,   311,   312,   312,   312,   313,   313,
1155           313,   314,   314,   314,   315,   315,   315,   315,   316,   316,
1156           316,   317,   317,   317,   318,   318,   318,   318,   319,   319,
1157           319,   320,   320,   320,   321,   321,   321,   321,   322,   322,
1158           322,   323,   323,   323,   324,   324,   324,   324,   325,   325,
1159           325,   326,   326,   326,   327,   327,   327,   327,   328,   328,
1160           328,   329,   329,   329,   330,   330,   330,   330,   331,   331,
1161           331,   332,   332,   332,   333,   333,   333,   333,   334,   334,
1162           334,   335,   335,   335,   336,   336,   336,   336,   337,   337,
1163           337,   338,   338,   338,   339,   339,   339,   339,   340,   340,
1164           340,   341,   341,   341,   342,   342,   342,   342,   343,   343,
1165           343,   344,   344,   344,   345,   345,   345,   346,   346,   346,
1166           346,   347,   347,   347,   348,   348,   348,   349,   349,   349,
1167           349,   350,   350,   350,   351,   351,   351,   352,   352,   352,
1168           352,   353,   353,   353,   354,   354,   354,   355,   355,   355,
1169           355,   356,   356,   356,   357,   357,   357,   358,   358,   358,
1170           358,   359,   359,   359,   360,   360,   360,   361,   361,   361,
1171           361,   362,   362,   362,   363,   363,   363,   364,   364,   364,
1172           364,   365,   365,   365,   366,   366,   366,   367,   367,   367,
1173           367,   368,   368,   368,   369,   369,   369,   370,   370,   370,
1174           370,   371,   371,   371,   372,   372,   372,   373,   373,   373,
1175           374,   374,   374,   374,   375,   375,   375,   376,   376,   376,
1176           377,   377,   377,   377,   378,   378,   378,   379,   379,   379,
1177           380,   380,   380,   380,   381,   381,   381,   382,   382,   382,
1178           383,   383,   383,   383,   384,   384,   384,   385,   385,   385,
1179           386,   386,   386,   386,   387,   387,   387,   388,   388,   388,
1180           389,   389,   389,   389,   390,   390,   390,   391,   391,   391,
1181           392,   392,   392,   392,   393,   393,   393,   394,   394,   394,
1182           395,   395,   395,   395,   396,   396,   396,   397,   397,   397,
1183           398,   398,   398,   398,   399,   399,   399,   400,   400,   400,
1184           401,   401,   401,   402,   402,   402,   402,   403,   403,   403,
1185           404,   404,   404,   405,   405,   405,   405,   406,   406,   406,
1186           407,   407,   407,   408,   408,   408,   408,   409,   409,   409,
1187           410,   410,   410,   411,   411,   411,   411,   412,   412,   412,
1188           413,   413,   413,   414,   414,   414,   414,   415,   415,   415,
1189           416,   416,   416,   417,   417,   417,   417,   418,   418,   418,
1190           419,   419,   419,   420,   420,   420,   420,   421,   421,   421,
1191           422,   422,   422,   423,   423,   423,   423,   424,   424,   424,
1192           425,   425,   425,   426,   426,   426,   426,   427,   427,   427,
1193           428,   428,   428,   429,   429,   429,   429,   430,   430,   430,
1194           431,   431,   431,   432,   432,   432,   433,   433,   433,   433,
1195           434,   434,   434,   435,   435,   435,   436,   436,   436,   436,
1196           437,   437,   437,   438,   438,   438,   439,   439,   439,   439,
1197           440,   440,   440,   441,   441,   441,   442,   442,   442,   442,
1198           443,   443,   443,   444,   444,   444,   445,   445,   445,   445,
1199           446,   446,   446,   447,   447,   447,   448,   448,   448,   448,
1200           449,   449,   449,   450,   450,   450,   451,   451,   451,   451,
1201           452,   452,   452,   453,   453,   453,   454,   454,   454,   454,
1202           455,   455,   455,   456,   456,   456,   457,   457,   457,   457,
1203           458,   458,   458,   459,   459,   459,   460,   460,   460,   461,
1204           461,   461,   461,   462,   462,   462,   463,   463,   463,   464,
1205           464,   464,   464,   465,   465,   465,   466,   466,   466,   467,
1206           467,   467,   467,   468,   468,   468,   469,   469,   469,   470,
1207           470,   470,   470,   471,   471,   471,   472,   472,   472,   473,
1208           473,   473,   473,   474,   474,   474,   475,   475,   475,   476,
1209           476,   476,   476,   477,   477,   477,   478,   478,   478,   479,
1210           479,   479,   479,   480,   480,   480,   481,   481,   481,   482,
1211           482,   482,   482,   483,   483,   483,   484,   484,   484,   485,
1212           485,   485,   485,   486,   486,   486,   487,   487,   487,   488,
1213           488,   488,   488,   489,   489,   489,   490,   490,   490,   491,
1214           491,   491,   492,   492,   492,   492,   493,   493,   493,   494,
1215           494,   494,   495,   495,   495,   495,   496,   496,   496,   497,
1216           497,   497,   498,   498,   498,   498,   499,   499,   499,   500,
1217           500,   500,   501,   501,   501,   501,   502,   502,   502,   503,
1218           503,   503,   504,   504,   504,   504,   505,   505,   505,   506,
1219           506,   506,   507,   507,   507,   507,   508,   508,   508,   509,
1220           509,   509,   510,   510,   510,   510,   511,   511,   511,   512,
1221           512,   512,   513,   513,   513,   513,   514,   514,   514,   515,
1222           515,   515,   516,   516,   516,   516,   517,   517,   517,   518,
1223           518,   518,   519,   519,   519,   520,   520,   520,   520,   521,
1224           521,   521,   522,   522,   522,   523,   523,   523,   523,   524,
1225           524,   524,   525,   525,   525,   526,   526,   526,   526,   527,
1226           527,   527,   528,   528,   528,   529,   529,   529,   529,   530,
1227           530,   530,   531,   531,   531,   532,   532,   532,   532,   533,
1228           533,   533,   534,   534,   534,   535,   535,   535,   535,   536,
1229           536,   536,   537,   537,   537,   538,   538,   538,   538,   539,
1230           539,   539,   540,   540,   540,   541,   541,   541,   541,   542,
1231           542,   542,   543,   543,   543,   544,   544,   544,   544,   545,
1232           545,   545,   546,   546,   546,   547,   547,   547,   548,   548,
1233           548,   548,   549,   549,   549,   550,   550,   550,   551,   551,
1234           551,   551,   552,   552,   552,   553,   553,   553,   554,   554,
1235           554,   554,   555,   555,   555,   556,   556,   556,   557,   557,
1236           557,   557,   558,   558,   558,   559,   559,   559,   560,   560,
1237           560,   560,   561,   561,   561,   562,   562,   562,   563,   563,
1238           563,   563,   564,   564,   564,   565,   565,   565,   566,   566,
1239           566,   566,   567,   567,   567,   568,   568,   568,   569,   569,
1240           569,   569,   570,   570,   570,   571,   571,   571,   572,   572,
1241           572,   572,   573,   573,   573,   574,   574,   574,   575,   575,
1242           575,   575,   576,   576,   576,   577,   577,   577,   578,   578,
1243           578,   579,   579,   579,   579,   580,   580,   580,   581,   581,
1244           581,   582,   582,   582,   582,   583,   583,   583,   584,   584,
1245           584,   585,   585,   585,   585,   586,   586,   586,   587,   587,
1246           587,   588,   588,   588,   588,   589,   589,   589,   590,   590,
1247           590,   591,   591,   591,   591,   592,   592,   592,   593,   593,
1248           593,   594,   594,   594,   594,   595,   595,   595,   596,   596,
1249           596,   597,   597,   597,   597,   598,   598,   598,   599,   599,
1250           599,   600,   600,   600,   600,   601,   601,   601,   602,   602,
1251           602,   603,   603,   603,   603,   604,   604,   604,   605,   605,
1252           605,   606,   606,   606,   607,   607,   607,   607,   608,   608,
1253           608,   609,   609,   609,   610,   610,   610,   610,   611,   611,
1254           611,   612,   612,   612,   613,   613,   613,   613,   614,   614,
1255           614,   615,   615,   615,   616,   616,   616,   616,   617,   617,
1256           617,   618,   618,   618,   619,   619,   619,   619,   620,   620,
1257           620,   621,   621,   621,   622,   622,   622,   622,   623,   623,
1258           623,   624,   624,   624,   625,   625,   625,   625,   626,   626,
1259           626,   627,   627,   627,   628,   628,   628,   628,   629,   629,
1260           629,   630,   630,   630,   631,   631,   631,   631,   632,   632,
1261           632,   633,   633,   633,   634,   634,   634,   634,   635,   635,
1262           635,   636,   636,   636,   637,   637,   637,   638,   638,   638,
1263           638,   639,   639,   639,   640,   640,   640,   641,   641,   641,
1264           641,   642,   642,   642,   643,   643,   643,   644,   644,   644,
1265           644,   645,   645,   645,   646,   646,   646,   647,   647,   647,
1266           647,   648,   648,   648,   649,   649,   649,   650,   650 };
1267  static ULLong pfive[27] = {
1268                 5ll,
1269                 25ll,
1270                 125ll,
1271                 625ll,
1272                 3125ll,
1273                 15625ll,
1274                 78125ll,
1275                 390625ll,
1276                 1953125ll,
1277                 9765625ll,
1278                 48828125ll,
1279                 244140625ll,
1280                 1220703125ll,
1281                 6103515625ll,
1282                 30517578125ll,
1283                 152587890625ll,
1284                 762939453125ll,
1285                 3814697265625ll,
1286                 19073486328125ll,
1287                 95367431640625ll,
1288                 476837158203125ll,
1289                 2384185791015625ll,
1290                 11920928955078125ll,
1291                 59604644775390625ll,
1292                 298023223876953125ll,
1293                 1490116119384765625ll,
1294                 7450580596923828125ll
1295                 };
1296
1297  static int pfivebits[25] = {3, 5, 7, 10, 12, 14, 17, 19, 21, 24, 26, 28, 31,
1298                              33, 35, 38, 40, 42, 45, 47, 49, 52, 54, 56, 59};
1299 #endif /*}*/
1300 #endif /*}} NO_LONG_LONG */
1301
1302 typedef union { double d; ULong L[2];
1303 #ifdef USE_BF96
1304         ULLong LL;
1305 #endif
1306         } U;
1307
1308 #ifdef IEEE_8087
1309 #define word0(x) (x)->L[1]
1310 #define word1(x) (x)->L[0]
1311 #else
1312 #define word0(x) (x)->L[0]
1313 #define word1(x) (x)->L[1]
1314 #endif
1315 #define dval(x) (x)->d
1316 #define LLval(x) (x)->LL
1317
1318 #ifndef STRTOD_DIGLIM
1319 #define STRTOD_DIGLIM 40
1320 #endif
1321
1322 #ifdef DIGLIM_DEBUG
1323 extern int strtod_diglim;
1324 #else
1325 #define strtod_diglim STRTOD_DIGLIM
1326 #endif
1327
1328 /* The following definition of Storeinc is appropriate for MIPS processors.
1329  * An alternative that might be better on some machines is
1330  * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
1331  */
1332 #if defined(IEEE_8087) + defined(VAX)
1333 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
1334 ((unsigned short *)a)[0] = (unsigned short)c, a++)
1335 #else
1336 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
1337 ((unsigned short *)a)[1] = (unsigned short)c, a++)
1338 #endif
1339
1340 /* #define P DBL_MANT_DIG */
1341 /* Ten_pmax = floor(P*log(2)/log(5)) */
1342 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
1343 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
1344 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
1345
1346 #ifdef IEEE_Arith
1347 #define Exp_shift  20
1348 #define Exp_shift1 20
1349 #define Exp_msk1    0x100000
1350 #define Exp_msk11   0x100000
1351 #define Exp_mask  0x7ff00000
1352 #define P 53
1353 #define Nbits 53
1354 #define Bias 1023
1355 #define Emax 1023
1356 #define Emin (-1022)
1357 #define Exp_1  0x3ff00000
1358 #define Exp_11 0x3ff00000
1359 #define Ebits 11
1360 #define Frac_mask  0xfffff
1361 #define Frac_mask1 0xfffff
1362 #define Ten_pmax 22
1363 #define Bletch 0x10
1364 #define Bndry_mask  0xfffff
1365 #define Bndry_mask1 0xfffff
1366 #define LSB 1
1367 #define Sign_bit 0x80000000
1368 #define Log2P 1
1369 #define Tiny0 0
1370 #define Tiny1 1
1371 #define Quick_max 14
1372 #define Int_max 14
1373 #ifndef NO_IEEE_Scale
1374 #define Avoid_Underflow
1375 #ifdef Flush_Denorm     /* debugging option */
1376 #undef Sudden_Underflow
1377 #endif
1378 #endif
1379
1380 #ifndef Flt_Rounds
1381 #ifdef FLT_ROUNDS
1382 #define Flt_Rounds FLT_ROUNDS
1383 #else
1384 #define Flt_Rounds 1
1385 #endif
1386 #endif /*Flt_Rounds*/
1387
1388 #ifdef Honor_FLT_ROUNDS
1389 #undef Check_FLT_ROUNDS
1390 #define Check_FLT_ROUNDS
1391 #else
1392 #define Rounding Flt_Rounds
1393 #endif
1394
1395 #else /* ifndef IEEE_Arith */
1396 #undef Check_FLT_ROUNDS
1397 #undef Honor_FLT_ROUNDS
1398 #undef SET_INEXACT
1399 #undef  Sudden_Underflow
1400 #define Sudden_Underflow
1401 #ifdef IBM
1402 #undef Flt_Rounds
1403 #define Flt_Rounds 0
1404 #define Exp_shift  24
1405 #define Exp_shift1 24
1406 #define Exp_msk1   0x1000000
1407 #define Exp_msk11  0x1000000
1408 #define Exp_mask  0x7f000000
1409 #define P 14
1410 #define Nbits 56
1411 #define Bias 65
1412 #define Emax 248
1413 #define Emin (-260)
1414 #define Exp_1  0x41000000
1415 #define Exp_11 0x41000000
1416 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
1417 #define Frac_mask  0xffffff
1418 #define Frac_mask1 0xffffff
1419 #define Bletch 4
1420 #define Ten_pmax 22
1421 #define Bndry_mask  0xefffff
1422 #define Bndry_mask1 0xffffff
1423 #define LSB 1
1424 #define Sign_bit 0x80000000
1425 #define Log2P 4
1426 #define Tiny0 0x100000
1427 #define Tiny1 0
1428 #define Quick_max 14
1429 #define Int_max 15
1430 #else /* VAX */
1431 #undef Flt_Rounds
1432 #define Flt_Rounds 1
1433 #define Exp_shift  23
1434 #define Exp_shift1 7
1435 #define Exp_msk1    0x80
1436 #define Exp_msk11   0x800000
1437 #define Exp_mask  0x7f80
1438 #define P 56
1439 #define Nbits 56
1440 #define Bias 129
1441 #define Emax 126
1442 #define Emin (-129)
1443 #define Exp_1  0x40800000
1444 #define Exp_11 0x4080
1445 #define Ebits 8
1446 #define Frac_mask  0x7fffff
1447 #define Frac_mask1 0xffff007f
1448 #define Ten_pmax 24
1449 #define Bletch 2
1450 #define Bndry_mask  0xffff007f
1451 #define Bndry_mask1 0xffff007f
1452 #define LSB 0x10000
1453 #define Sign_bit 0x8000
1454 #define Log2P 1
1455 #define Tiny0 0x80
1456 #define Tiny1 0
1457 #define Quick_max 15
1458 #define Int_max 15
1459 #endif /* IBM, VAX */
1460 #endif /* IEEE_Arith */
1461
1462 #ifndef IEEE_Arith
1463 #define ROUND_BIASED
1464 #else
1465 #ifdef ROUND_BIASED_without_Round_Up
1466 #undef  ROUND_BIASED
1467 #define ROUND_BIASED
1468 #endif
1469 #endif
1470
1471 #ifdef RND_PRODQUOT
1472 #define rounded_product(a,b) a = rnd_prod(a, b)
1473 #define rounded_quotient(a,b) a = rnd_quot(a, b)
1474 extern double rnd_prod(double, double), rnd_quot(double, double);
1475 #else
1476 #define rounded_product(a,b) a *= b
1477 #define rounded_quotient(a,b) a /= b
1478 #endif
1479
1480 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
1481 #define Big1 0xffffffff
1482
1483 #ifndef Pack_32
1484 #define Pack_32
1485 #endif
1486
1487 typedef struct BCinfo BCinfo;
1488  struct
1489 BCinfo { int dp0, dp1, dplen, dsign, e0, inexact, nd, nd0, rounding, scale, uflchk; };
1490
1491 #define FFFFFFFF 0xffffffffUL
1492
1493 #ifdef MULTIPLE_THREADS
1494 #define MTa , PTI
1495 #define MTb , &TI
1496 #define MTd , ThInfo **PTI
1497 static unsigned int maxthreads = 0;
1498 #else
1499 #define MTa /*nothing*/
1500 #define MTb /*nothing*/
1501 #define MTd /*nothing*/
1502 #endif
1503
1504 #define Kmax 7
1505
1506 #ifdef __cplusplus
1507 extern "C" double strtod(const char *s00, char **se);
1508 extern "C" char *dtoa(double d, int mode, int ndigits,
1509                         int *decpt, int *sign, char **rve);
1510 #endif
1511
1512  struct
1513 Bigint {
1514         struct Bigint *next;
1515         int k, maxwds, sign, wds;
1516         ULong x[1];
1517         };
1518
1519  typedef struct Bigint Bigint;
1520  typedef struct
1521 ThInfo {
1522         Bigint *Freelist[Kmax+1];
1523         Bigint *P5s;
1524         } ThInfo;
1525
1526  static ThInfo TI0;
1527
1528 #ifdef MULTIPLE_THREADS
1529  static ThInfo *TI1;
1530  static int TI0_used;
1531
1532  void
1533 set_max_dtoa_threads(unsigned int n)
1534 {
1535         size_t L;
1536
1537         if (n > maxthreads) {
1538                 L = n*sizeof(ThInfo);
1539                 if (TI1) {
1540                         TI1 = (ThInfo*)REALLOC(TI1, L);
1541                         memset(TI1 + maxthreads, 0, (n-maxthreads)*sizeof(ThInfo));
1542                         }
1543                 else {
1544                         TI1 = (ThInfo*)MALLOC(L);
1545                         if (TI0_used) {
1546                                 memcpy(TI1, &TI0, sizeof(ThInfo));
1547                                 if (n > 1)
1548                                         memset(TI1 + 1, 0, L - sizeof(ThInfo));
1549                                 memset(&TI0, 0, sizeof(ThInfo));
1550                                 }
1551                         else
1552                                 memset(TI1, 0, L);
1553                         }
1554                 maxthreads = n;
1555                 }
1556         }
1557
1558  static ThInfo*
1559 get_TI(void)
1560 {
1561         unsigned int thno = dtoa_get_threadno();
1562         if (thno < maxthreads)
1563                 return TI1 + thno;
1564         if (thno == 0)
1565                 TI0_used = 1;
1566         return &TI0;
1567         }
1568 #define freelist TI->Freelist
1569 #define p5s TI->P5s
1570 #else
1571 #define freelist TI0.Freelist
1572 #define p5s TI0.P5s
1573 #endif
1574
1575  static Bigint *
1576 Balloc(int k MTd)
1577 {
1578         int x;
1579         Bigint *rv;
1580 #ifndef Omit_Private_Memory
1581         unsigned int len;
1582 #endif
1583 #ifdef MULTIPLE_THREADS
1584         ThInfo *TI;
1585
1586         if (!(TI = *PTI))
1587                 *PTI = TI = get_TI();
1588         if (TI == &TI0)
1589                 ACQUIRE_DTOA_LOCK(0);
1590 #endif
1591         /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
1592         /* but this case seems very unlikely. */
1593         if (k <= Kmax && (rv = freelist[k]))
1594                 freelist[k] = rv->next;
1595         else {
1596                 x = 1 << k;
1597 #ifdef Omit_Private_Memory
1598                 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
1599 #else
1600                 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
1601                         /sizeof(double);
1602                 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem
1603 #ifdef MULTIPLE_THREADS
1604                         && TI == TI1
1605 #endif
1606                         ) {
1607                         rv = (Bigint*)pmem_next;
1608                         pmem_next += len;
1609                         }
1610                 else
1611                         rv = (Bigint*)MALLOC(len*sizeof(double));
1612 #endif
1613                 rv->k = k;
1614                 rv->maxwds = x;
1615                 }
1616 #ifdef MULTIPLE_THREADS
1617         if (TI == &TI0)
1618                 FREE_DTOA_LOCK(0);
1619 #endif
1620         rv->sign = rv->wds = 0;
1621         return rv;
1622         }
1623
1624  static void
1625 Bfree(Bigint *v MTd)
1626 {
1627 #ifdef MULTIPLE_THREADS
1628         ThInfo *TI;
1629 #endif
1630         if (v) {
1631                 if (v->k > Kmax)
1632                         FREE((void*)v);
1633                 else {
1634 #ifdef MULTIPLE_THREADS
1635                         if (!(TI = *PTI))
1636                                 *PTI = TI = get_TI();
1637                         if (TI == &TI0)
1638                                 ACQUIRE_DTOA_LOCK(0);
1639 #endif
1640                         v->next = freelist[v->k];
1641                         freelist[v->k] = v;
1642 #ifdef MULTIPLE_THREADS
1643                         if (TI == &TI0)
1644                                 FREE_DTOA_LOCK(0);
1645 #endif
1646                         }
1647                 }
1648         }
1649
1650 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
1651 y->wds*sizeof(Long) + 2*sizeof(int))
1652
1653  static Bigint *
1654 multadd(Bigint *b, int m, int a MTd)    /* multiply by m and add a */
1655 {
1656         int i, wds;
1657 #ifdef ULLong
1658         ULong *x;
1659         ULLong carry, y;
1660 #else
1661         ULong carry, *x, y;
1662 #ifdef Pack_32
1663         ULong xi, z;
1664 #endif
1665 #endif
1666         Bigint *b1;
1667
1668         wds = b->wds;
1669         x = b->x;
1670         i = 0;
1671         carry = a;
1672         do {
1673 #ifdef ULLong
1674                 y = *x * (ULLong)m + carry;
1675                 carry = y >> 32;
1676                 *x++ = y & FFFFFFFF;
1677 #else
1678 #ifdef Pack_32
1679                 xi = *x;
1680                 y = (xi & 0xffff) * m + carry;
1681                 z = (xi >> 16) * m + (y >> 16);
1682                 carry = z >> 16;
1683                 *x++ = (z << 16) + (y & 0xffff);
1684 #else
1685                 y = *x * m + carry;
1686                 carry = y >> 16;
1687                 *x++ = y & 0xffff;
1688 #endif
1689 #endif
1690                 }
1691                 while(++i < wds);
1692         if (carry) {
1693                 if (wds >= b->maxwds) {
1694                         b1 = Balloc(b->k+1 MTa);
1695                         Bcopy(b1, b);
1696                         Bfree(b MTa);
1697                         b = b1;
1698                         }
1699                 b->x[wds++] = carry;
1700                 b->wds = wds;
1701                 }
1702         return b;
1703         }
1704
1705  static Bigint *
1706 s2b(const char *s, int nd0, int nd, ULong y9, int dplen MTd)
1707 {
1708         Bigint *b;
1709         int i, k;
1710         Long x, y;
1711
1712         x = (nd + 8) / 9;
1713         for(k = 0, y = 1; x > y; y <<= 1, k++) ;
1714 #ifdef Pack_32
1715         b = Balloc(k MTa);
1716         b->x[0] = y9;
1717         b->wds = 1;
1718 #else
1719         b = Balloc(k+1 MTa);
1720         b->x[0] = y9 & 0xffff;
1721         b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
1722 #endif
1723
1724         i = 9;
1725         if (9 < nd0) {
1726                 s += 9;
1727                 do b = multadd(b, 10, *s++ - '0' MTa);
1728                         while(++i < nd0);
1729                 s += dplen;
1730                 }
1731         else
1732                 s += dplen + 9;
1733         for(; i < nd; i++)
1734                 b = multadd(b, 10, *s++ - '0' MTa);
1735         return b;
1736         }
1737
1738  static int
1739 hi0bits(ULong x)
1740 {
1741         int k = 0;
1742
1743         if (!(x & 0xffff0000)) {
1744                 k = 16;
1745                 x <<= 16;
1746                 }
1747         if (!(x & 0xff000000)) {
1748                 k += 8;
1749                 x <<= 8;
1750                 }
1751         if (!(x & 0xf0000000)) {
1752                 k += 4;
1753                 x <<= 4;
1754                 }
1755         if (!(x & 0xc0000000)) {
1756                 k += 2;
1757                 x <<= 2;
1758                 }
1759         if (!(x & 0x80000000)) {
1760                 k++;
1761                 if (!(x & 0x40000000))
1762                         return 32;
1763                 }
1764         return k;
1765         }
1766
1767  static int
1768 lo0bits(ULong *y)
1769 {
1770         int k;
1771         ULong x = *y;
1772
1773         if (x & 7) {
1774                 if (x & 1)
1775                         return 0;
1776                 if (x & 2) {
1777                         *y = x >> 1;
1778                         return 1;
1779                         }
1780                 *y = x >> 2;
1781                 return 2;
1782                 }
1783         k = 0;
1784         if (!(x & 0xffff)) {
1785                 k = 16;
1786                 x >>= 16;
1787                 }
1788         if (!(x & 0xff)) {
1789                 k += 8;
1790                 x >>= 8;
1791                 }
1792         if (!(x & 0xf)) {
1793                 k += 4;
1794                 x >>= 4;
1795                 }
1796         if (!(x & 0x3)) {
1797                 k += 2;
1798                 x >>= 2;
1799                 }
1800         if (!(x & 1)) {
1801                 k++;
1802                 x >>= 1;
1803                 if (!x)
1804                         return 32;
1805                 }
1806         *y = x;
1807         return k;
1808         }
1809
1810  static Bigint *
1811 i2b(int i MTd)
1812 {
1813         Bigint *b;
1814
1815         b = Balloc(1 MTa);
1816         b->x[0] = i;
1817         b->wds = 1;
1818         return b;
1819         }
1820
1821  static Bigint *
1822 mult(Bigint *a, Bigint *b MTd)
1823 {
1824         Bigint *c;
1825         int k, wa, wb, wc;
1826         ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
1827         ULong y;
1828 #ifdef ULLong
1829         ULLong carry, z;
1830 #else
1831         ULong carry, z;
1832 #ifdef Pack_32
1833         ULong z2;
1834 #endif
1835 #endif
1836
1837         if (a->wds < b->wds) {
1838                 c = a;
1839                 a = b;
1840                 b = c;
1841                 }
1842         k = a->k;
1843         wa = a->wds;
1844         wb = b->wds;
1845         wc = wa + wb;
1846         if (wc > a->maxwds)
1847                 k++;
1848         c = Balloc(k MTa);
1849         for(x = c->x, xa = x + wc; x < xa; x++)
1850                 *x = 0;
1851         xa = a->x;
1852         xae = xa + wa;
1853         xb = b->x;
1854         xbe = xb + wb;
1855         xc0 = c->x;
1856 #ifdef ULLong
1857         for(; xb < xbe; xc0++) {
1858                 if ((y = *xb++)) {
1859                         x = xa;
1860                         xc = xc0;
1861                         carry = 0;
1862                         do {
1863                                 z = *x++ * (ULLong)y + *xc + carry;
1864                                 carry = z >> 32;
1865                                 *xc++ = z & FFFFFFFF;
1866                                 }
1867                                 while(x < xae);
1868                         *xc = carry;
1869                         }
1870                 }
1871 #else
1872 #ifdef Pack_32
1873         for(; xb < xbe; xb++, xc0++) {
1874                 if (y = *xb & 0xffff) {
1875                         x = xa;
1876                         xc = xc0;
1877                         carry = 0;
1878                         do {
1879                                 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
1880                                 carry = z >> 16;
1881                                 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
1882                                 carry = z2 >> 16;
1883                                 Storeinc(xc, z2, z);
1884                                 }
1885                                 while(x < xae);
1886                         *xc = carry;
1887                         }
1888                 if (y = *xb >> 16) {
1889                         x = xa;
1890                         xc = xc0;
1891                         carry = 0;
1892                         z2 = *xc;
1893                         do {
1894                                 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
1895                                 carry = z >> 16;
1896                                 Storeinc(xc, z, z2);
1897                                 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
1898                                 carry = z2 >> 16;
1899                                 }
1900                                 while(x < xae);
1901                         *xc = z2;
1902                         }
1903                 }
1904 #else
1905         for(; xb < xbe; xc0++) {
1906                 if (y = *xb++) {
1907                         x = xa;
1908                         xc = xc0;
1909                         carry = 0;
1910                         do {
1911                                 z = *x++ * y + *xc + carry;
1912                                 carry = z >> 16;
1913                                 *xc++ = z & 0xffff;
1914                                 }
1915                                 while(x < xae);
1916                         *xc = carry;
1917                         }
1918                 }
1919 #endif
1920 #endif
1921         for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
1922         c->wds = wc;
1923         return c;
1924         }
1925
1926  static Bigint *
1927 pow5mult(Bigint *b, int k MTd)
1928 {
1929         Bigint *b1, *p5, *p51;
1930 #ifdef MULTIPLE_THREADS
1931         ThInfo *TI;
1932 #endif
1933         int i;
1934         static int p05[3] = { 5, 25, 125 };
1935
1936         if ((i = k & 3))
1937                 b = multadd(b, p05[i-1], 0 MTa);
1938
1939         if (!(k >>= 2))
1940                 return b;
1941 #ifdef  MULTIPLE_THREADS
1942         if (!(TI = *PTI))
1943                 *PTI = TI = get_TI();
1944 #endif
1945         if (!(p5 = p5s)) {
1946                 /* first time */
1947 #ifdef MULTIPLE_THREADS
1948                 if (!(TI = *PTI))
1949                         *PTI = TI = get_TI();
1950                 if (TI == &TI0)
1951                         ACQUIRE_DTOA_LOCK(1);
1952                 if (!(p5 = p5s)) {
1953                         p5 = p5s = i2b(625 MTa);
1954                         p5->next = 0;
1955                         }
1956                 if (TI == &TI0)
1957                         FREE_DTOA_LOCK(1);
1958 #else
1959                 p5 = p5s = i2b(625 MTa);
1960                 p5->next = 0;
1961 #endif
1962                 }
1963         for(;;) {
1964                 if (k & 1) {
1965                         b1 = mult(b, p5 MTa);
1966                         Bfree(b MTa);
1967                         b = b1;
1968                         }
1969                 if (!(k >>= 1))
1970                         break;
1971                 if (!(p51 = p5->next)) {
1972 #ifdef MULTIPLE_THREADS
1973                         if (!TI && !(TI = *PTI))
1974                                 *PTI = TI = get_TI();
1975                         if (TI == &TI0)
1976                                 ACQUIRE_DTOA_LOCK(1);
1977                         if (!(p51 = p5->next)) {
1978                                 p51 = p5->next = mult(p5,p5 MTa);
1979                                 p51->next = 0;
1980                                 }
1981                         if (TI == &TI0)
1982                                 FREE_DTOA_LOCK(1);
1983 #else
1984                         p51 = p5->next = mult(p5,p5);
1985                         p51->next = 0;
1986 #endif
1987                         }
1988                 p5 = p51;
1989                 }
1990         return b;
1991         }
1992
1993  static Bigint *
1994 lshift(Bigint *b, int k MTd)
1995 {
1996         int i, k1, n, n1;
1997         Bigint *b1;
1998         ULong *x, *x1, *xe, z;
1999
2000 #ifdef Pack_32
2001         n = k >> 5;
2002 #else
2003         n = k >> 4;
2004 #endif
2005         k1 = b->k;
2006         n1 = n + b->wds + 1;
2007         for(i = b->maxwds; n1 > i; i <<= 1)
2008                 k1++;
2009         b1 = Balloc(k1 MTa);
2010         x1 = b1->x;
2011         for(i = 0; i < n; i++)
2012                 *x1++ = 0;
2013         x = b->x;
2014         xe = x + b->wds;
2015 #ifdef Pack_32
2016         if (k &= 0x1f) {
2017                 k1 = 32 - k;
2018                 z = 0;
2019                 do {
2020                         *x1++ = *x << k | z;
2021                         z = *x++ >> k1;
2022                         }
2023                         while(x < xe);
2024                 if ((*x1 = z))
2025                         ++n1;
2026                 }
2027 #else
2028         if (k &= 0xf) {
2029                 k1 = 16 - k;
2030                 z = 0;
2031                 do {
2032                         *x1++ = *x << k  & 0xffff | z;
2033                         z = *x++ >> k1;
2034                         }
2035                         while(x < xe);
2036                 if (*x1 = z)
2037                         ++n1;
2038                 }
2039 #endif
2040         else do
2041                 *x1++ = *x++;
2042                 while(x < xe);
2043         b1->wds = n1 - 1;
2044         Bfree(b MTa);
2045         return b1;
2046         }
2047
2048  static int
2049 cmp(Bigint *a, Bigint *b)
2050 {
2051         ULong *xa, *xa0, *xb, *xb0;
2052         int i, j;
2053
2054         i = a->wds;
2055         j = b->wds;
2056 #ifdef DEBUG
2057         if (i > 1 && !a->x[i-1])
2058                 Bug("cmp called with a->x[a->wds-1] == 0");
2059         if (j > 1 && !b->x[j-1])
2060                 Bug("cmp called with b->x[b->wds-1] == 0");
2061 #endif
2062         if (i -= j)
2063                 return i;
2064         xa0 = a->x;
2065         xa = xa0 + j;
2066         xb0 = b->x;
2067         xb = xb0 + j;
2068         for(;;) {
2069                 if (*--xa != *--xb)
2070                         return *xa < *xb ? -1 : 1;
2071                 if (xa <= xa0)
2072                         break;
2073                 }
2074         return 0;
2075         }
2076
2077  static Bigint *
2078 diff(Bigint *a, Bigint *b MTd)
2079 {
2080         Bigint *c;
2081         int i, wa, wb;
2082         ULong *xa, *xae, *xb, *xbe, *xc;
2083 #ifdef ULLong
2084         ULLong borrow, y;
2085 #else
2086         ULong borrow, y;
2087 #ifdef Pack_32
2088         ULong z;
2089 #endif
2090 #endif
2091
2092         i = cmp(a,b);
2093         if (!i) {
2094                 c = Balloc(0 MTa);
2095                 c->wds = 1;
2096                 c->x[0] = 0;
2097                 return c;
2098                 }
2099         if (i < 0) {
2100                 c = a;
2101                 a = b;
2102                 b = c;
2103                 i = 1;
2104                 }
2105         else
2106                 i = 0;
2107         c = Balloc(a->k MTa);
2108         c->sign = i;
2109         wa = a->wds;
2110         xa = a->x;
2111         xae = xa + wa;
2112         wb = b->wds;
2113         xb = b->x;
2114         xbe = xb + wb;
2115         xc = c->x;
2116         borrow = 0;
2117 #ifdef ULLong
2118         do {
2119                 y = (ULLong)*xa++ - *xb++ - borrow;
2120                 borrow = y >> 32 & (ULong)1;
2121                 *xc++ = y & FFFFFFFF;
2122                 }
2123                 while(xb < xbe);
2124         while(xa < xae) {
2125                 y = *xa++ - borrow;
2126                 borrow = y >> 32 & (ULong)1;
2127                 *xc++ = y & FFFFFFFF;
2128                 }
2129 #else
2130 #ifdef Pack_32
2131         do {
2132                 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
2133                 borrow = (y & 0x10000) >> 16;
2134                 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
2135                 borrow = (z & 0x10000) >> 16;
2136                 Storeinc(xc, z, y);
2137                 }
2138                 while(xb < xbe);
2139         while(xa < xae) {
2140                 y = (*xa & 0xffff) - borrow;
2141                 borrow = (y & 0x10000) >> 16;
2142                 z = (*xa++ >> 16) - borrow;
2143                 borrow = (z & 0x10000) >> 16;
2144                 Storeinc(xc, z, y);
2145                 }
2146 #else
2147         do {
2148                 y = *xa++ - *xb++ - borrow;
2149                 borrow = (y & 0x10000) >> 16;
2150                 *xc++ = y & 0xffff;
2151                 }
2152                 while(xb < xbe);
2153         while(xa < xae) {
2154                 y = *xa++ - borrow;
2155                 borrow = (y & 0x10000) >> 16;
2156                 *xc++ = y & 0xffff;
2157                 }
2158 #endif
2159 #endif
2160         while(!*--xc)
2161                 wa--;
2162         c->wds = wa;
2163         return c;
2164         }
2165
2166  static double
2167 ulp(U *x)
2168 {
2169         Long L;
2170         U u;
2171
2172         L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
2173 #ifndef Avoid_Underflow
2174 #ifndef Sudden_Underflow
2175         if (L > 0) {
2176 #endif
2177 #endif
2178 #ifdef IBM
2179                 L |= Exp_msk1 >> 4;
2180 #endif
2181                 word0(&u) = L;
2182                 word1(&u) = 0;
2183 #ifndef Avoid_Underflow
2184 #ifndef Sudden_Underflow
2185                 }
2186         else {
2187                 L = -L >> Exp_shift;
2188                 if (L < Exp_shift) {
2189                         word0(&u) = 0x80000 >> L;
2190                         word1(&u) = 0;
2191                         }
2192                 else {
2193                         word0(&u) = 0;
2194                         L -= Exp_shift;
2195                         word1(&u) = L >= 31 ? 1 : 1 << 31 - L;
2196                         }
2197                 }
2198 #endif
2199 #endif
2200         return dval(&u);
2201         }
2202
2203  static double
2204 b2d(Bigint *a, int *e)
2205 {
2206         ULong *xa, *xa0, w, y, z;
2207         int k;
2208         U d;
2209 #ifdef VAX
2210         ULong d0, d1;
2211 #else
2212 #define d0 word0(&d)
2213 #define d1 word1(&d)
2214 #endif
2215
2216         xa0 = a->x;
2217         xa = xa0 + a->wds;
2218         y = *--xa;
2219 #ifdef DEBUG
2220         if (!y) Bug("zero y in b2d");
2221 #endif
2222         k = hi0bits(y);
2223         *e = 32 - k;
2224 #ifdef Pack_32
2225         if (k < Ebits) {
2226                 d0 = Exp_1 | y >> (Ebits - k);
2227                 w = xa > xa0 ? *--xa : 0;
2228                 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
2229                 goto ret_d;
2230                 }
2231         z = xa > xa0 ? *--xa : 0;
2232         if (k -= Ebits) {
2233                 d0 = Exp_1 | y << k | z >> (32 - k);
2234                 y = xa > xa0 ? *--xa : 0;
2235                 d1 = z << k | y >> (32 - k);
2236                 }
2237         else {
2238                 d0 = Exp_1 | y;
2239                 d1 = z;
2240                 }
2241 #else
2242         if (k < Ebits + 16) {
2243                 z = xa > xa0 ? *--xa : 0;
2244                 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
2245                 w = xa > xa0 ? *--xa : 0;
2246                 y = xa > xa0 ? *--xa : 0;
2247                 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
2248                 goto ret_d;
2249                 }
2250         z = xa > xa0 ? *--xa : 0;
2251         w = xa > xa0 ? *--xa : 0;
2252         k -= Ebits + 16;
2253         d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
2254         y = xa > xa0 ? *--xa : 0;
2255         d1 = w << k + 16 | y << k;
2256 #endif
2257  ret_d:
2258 #ifdef VAX
2259         word0(&d) = d0 >> 16 | d0 << 16;
2260         word1(&d) = d1 >> 16 | d1 << 16;
2261 #else
2262 #undef d0
2263 #undef d1
2264 #endif
2265         return dval(&d);
2266         }
2267
2268  static Bigint *
2269 d2b(U *d, int *e, int *bits MTd)
2270 {
2271         Bigint *b;
2272         int de, k;
2273         ULong *x, y, z;
2274 #ifndef Sudden_Underflow
2275         int i;
2276 #endif
2277 #ifdef VAX
2278         ULong d0, d1;
2279         d0 = word0(d) >> 16 | word0(d) << 16;
2280         d1 = word1(d) >> 16 | word1(d) << 16;
2281 #else
2282 #define d0 word0(d)
2283 #define d1 word1(d)
2284 #endif
2285
2286 #ifdef Pack_32
2287         b = Balloc(1 MTa);
2288 #else
2289         b = Balloc(2 MTa);
2290 #endif
2291         x = b->x;
2292
2293         z = d0 & Frac_mask;
2294         d0 &= 0x7fffffff;       /* clear sign bit, which we ignore */
2295 #ifdef Sudden_Underflow
2296         de = (int)(d0 >> Exp_shift);
2297 #ifndef IBM
2298         z |= Exp_msk11;
2299 #endif
2300 #else
2301         if ((de = (int)(d0 >> Exp_shift)))
2302                 z |= Exp_msk1;
2303 #endif
2304 #ifdef Pack_32
2305         if ((y = d1)) {
2306                 if ((k = lo0bits(&y))) {
2307                         x[0] = y | z << (32 - k);
2308                         z >>= k;
2309                         }
2310                 else
2311                         x[0] = y;
2312 #ifndef Sudden_Underflow
2313                 i =
2314 #endif
2315                     b->wds = (x[1] = z) ? 2 : 1;
2316                 }
2317         else {
2318                 k = lo0bits(&z);
2319                 x[0] = z;
2320 #ifndef Sudden_Underflow
2321                 i =
2322 #endif
2323                     b->wds = 1;
2324                 k += 32;
2325                 }
2326 #else
2327         if (y = d1) {
2328                 if (k = lo0bits(&y))
2329                         if (k >= 16) {
2330                                 x[0] = y | z << 32 - k & 0xffff;
2331                                 x[1] = z >> k - 16 & 0xffff;
2332                                 x[2] = z >> k;
2333                                 i = 2;
2334                                 }
2335                         else {
2336                                 x[0] = y & 0xffff;
2337                                 x[1] = y >> 16 | z << 16 - k & 0xffff;
2338                                 x[2] = z >> k & 0xffff;
2339                                 x[3] = z >> k+16;
2340                                 i = 3;
2341                                 }
2342                 else {
2343                         x[0] = y & 0xffff;
2344                         x[1] = y >> 16;
2345                         x[2] = z & 0xffff;
2346                         x[3] = z >> 16;
2347                         i = 3;
2348                         }
2349                 }
2350         else {
2351 #ifdef DEBUG
2352                 if (!z)
2353                         Bug("Zero passed to d2b");
2354 #endif
2355                 k = lo0bits(&z);
2356                 if (k >= 16) {
2357                         x[0] = z;
2358                         i = 0;
2359                         }
2360                 else {
2361                         x[0] = z & 0xffff;
2362                         x[1] = z >> 16;
2363                         i = 1;
2364                         }
2365                 k += 32;
2366                 }
2367         while(!x[i])
2368                 --i;
2369         b->wds = i + 1;
2370 #endif
2371 #ifndef Sudden_Underflow
2372         if (de) {
2373 #endif
2374 #ifdef IBM
2375                 *e = (de - Bias - (P-1) << 2) + k;
2376                 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
2377 #else
2378                 *e = de - Bias - (P-1) + k;
2379                 *bits = P - k;
2380 #endif
2381 #ifndef Sudden_Underflow
2382                 }
2383         else {
2384                 *e = de - Bias - (P-1) + 1 + k;
2385 #ifdef Pack_32
2386                 *bits = 32*i - hi0bits(x[i-1]);
2387 #else
2388                 *bits = (i+2)*16 - hi0bits(x[i]);
2389 #endif
2390                 }
2391 #endif
2392         return b;
2393         }
2394 #undef d0
2395 #undef d1
2396
2397  static double
2398 ratio(Bigint *a, Bigint *b)
2399 {
2400         U da, db;
2401         int k, ka, kb;
2402
2403         dval(&da) = b2d(a, &ka);
2404         dval(&db) = b2d(b, &kb);
2405 #ifdef Pack_32
2406         k = ka - kb + 32*(a->wds - b->wds);
2407 #else
2408         k = ka - kb + 16*(a->wds - b->wds);
2409 #endif
2410 #ifdef IBM
2411         if (k > 0) {
2412                 word0(&da) += (k >> 2)*Exp_msk1;
2413                 if (k &= 3)
2414                         dval(&da) *= 1 << k;
2415                 }
2416         else {
2417                 k = -k;
2418                 word0(&db) += (k >> 2)*Exp_msk1;
2419                 if (k &= 3)
2420                         dval(&db) *= 1 << k;
2421                 }
2422 #else
2423         if (k > 0)
2424                 word0(&da) += k*Exp_msk1;
2425         else {
2426                 k = -k;
2427                 word0(&db) += k*Exp_msk1;
2428                 }
2429 #endif
2430         return dval(&da) / dval(&db);
2431         }
2432
2433  static const double
2434 tens[] = {
2435                 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
2436                 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
2437                 1e20, 1e21, 1e22
2438 #ifdef VAX
2439                 , 1e23, 1e24
2440 #endif
2441                 };
2442
2443  static const double
2444 #ifdef IEEE_Arith
2445 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
2446 static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
2447 #ifdef Avoid_Underflow
2448                 9007199254740992.*9007199254740992.e-256
2449                 /* = 2^106 * 1e-256 */
2450 #else
2451                 1e-256
2452 #endif
2453                 };
2454 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
2455 /* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
2456 #define Scale_Bit 0x10
2457 #define n_bigtens 5
2458 #else
2459 #ifdef IBM
2460 bigtens[] = { 1e16, 1e32, 1e64 };
2461 static const double tinytens[] = { 1e-16, 1e-32, 1e-64 };
2462 #define n_bigtens 3
2463 #else
2464 bigtens[] = { 1e16, 1e32 };
2465 static const double tinytens[] = { 1e-16, 1e-32 };
2466 #define n_bigtens 2
2467 #endif
2468 #endif
2469
2470 #undef Need_Hexdig
2471 #ifdef INFNAN_CHECK
2472 #ifndef No_Hex_NaN
2473 #define Need_Hexdig
2474 #endif
2475 #endif
2476
2477 #ifndef Need_Hexdig
2478 #ifndef NO_HEX_FP
2479 #define Need_Hexdig
2480 #endif
2481 #endif
2482
2483 #ifdef Need_Hexdig /*{*/
2484 #if 0
2485 static unsigned char hexdig[256];
2486
2487  static void
2488 htinit(unsigned char *h, unsigned char *s, int inc)
2489 {
2490         int i, j;
2491         for(i = 0; (j = s[i]) !=0; i++)
2492                 h[j] = i + inc;
2493         }
2494
2495  static void
2496 hexdig_init(void)       /* Use of hexdig_init omitted 20121220 to avoid a */
2497                         /* race condition when multiple threads are used. */
2498 {
2499 #define USC (unsigned char *)
2500         htinit(hexdig, USC "0123456789", 0x10);
2501         htinit(hexdig, USC "abcdef", 0x10 + 10);
2502         htinit(hexdig, USC "ABCDEF", 0x10 + 10);
2503         }
2504 #else
2505 static unsigned char hexdig[256] = {
2506         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2507         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2508         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2509         16,17,18,19,20,21,22,23,24,25,0,0,0,0,0,0,
2510         0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0,
2511         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2512         0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0,
2513         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2514         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2515         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2516         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2517         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2518         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2519         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2520         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2521         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
2522         };
2523 #endif
2524 #endif /* } Need_Hexdig */
2525
2526 #ifdef INFNAN_CHECK
2527
2528 #ifndef NAN_WORD0
2529 #define NAN_WORD0 0x7ff80000
2530 #endif
2531
2532 #ifndef NAN_WORD1
2533 #define NAN_WORD1 0
2534 #endif
2535
2536  static int
2537 match(const char **sp, const char *t)
2538 {
2539         int c, d;
2540         const char *s = *sp;
2541
2542         while((d = *t++)) {
2543                 if ((c = *++s) >= 'A' && c <= 'Z')
2544                         c += 'a' - 'A';
2545                 if (c != d)
2546                         return 0;
2547                 }
2548         *sp = s + 1;
2549         return 1;
2550         }
2551
2552 #ifndef No_Hex_NaN
2553  static void
2554 hexnan(U *rvp, const char **sp)
2555 {
2556         ULong c, x[2];
2557         const char *s;
2558         int c1, havedig, udx0, xshift;
2559
2560         /**** if (!hexdig['0']) hexdig_init(); ****/
2561         x[0] = x[1] = 0;
2562         havedig = xshift = 0;
2563         udx0 = 1;
2564         s = *sp;
2565         /* allow optional initial 0x or 0X */
2566         while((c = *(const unsigned char*)(s+1)) && c <= ' ')
2567                 ++s;
2568         if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X'))
2569                 s += 2;
2570         while((c = *(const unsigned char*)++s)) {
2571                 if ((c1 = hexdig[c]))
2572                         c  = c1 & 0xf;
2573                 else if (c <= ' ') {
2574                         if (udx0 && havedig) {
2575                                 udx0 = 0;
2576                                 xshift = 1;
2577                                 }
2578                         continue;
2579                         }
2580 #ifdef GDTOA_NON_PEDANTIC_NANCHECK
2581                 else if (/*(*/ c == ')' && havedig) {
2582                         *sp = s + 1;
2583                         break;
2584                         }
2585                 else
2586                         return; /* invalid form: don't change *sp */
2587 #else
2588                 else {
2589                         do {
2590                                 if (/*(*/ c == ')') {
2591                                         *sp = s + 1;
2592                                         break;
2593                                         }
2594                                 } while((c = *++s));
2595                         break;
2596                         }
2597 #endif
2598                 havedig = 1;
2599                 if (xshift) {
2600                         xshift = 0;
2601                         x[0] = x[1];
2602                         x[1] = 0;
2603                         }
2604                 if (udx0)
2605                         x[0] = (x[0] << 4) | (x[1] >> 28);
2606                 x[1] = (x[1] << 4) | c;
2607                 }
2608         if ((x[0] &= 0xfffff) || x[1]) {
2609                 word0(rvp) = Exp_mask | x[0];
2610                 word1(rvp) = x[1];
2611                 }
2612         }
2613 #endif /*No_Hex_NaN*/
2614 #endif /* INFNAN_CHECK */
2615
2616 #ifdef Pack_32
2617 #define ULbits 32
2618 #define kshift 5
2619 #define kmask 31
2620 #else
2621 #define ULbits 16
2622 #define kshift 4
2623 #define kmask 15
2624 #endif
2625
2626 #if !defined(NO_HEX_FP) || defined(Honor_FLT_ROUNDS) /*{*/
2627  static Bigint *
2628 increment(Bigint *b MTd)
2629 {
2630         ULong *x, *xe;
2631         Bigint *b1;
2632
2633         x = b->x;
2634         xe = x + b->wds;
2635         do {
2636                 if (*x < (ULong)0xffffffffL) {
2637                         ++*x;
2638                         return b;
2639                         }
2640                 *x++ = 0;
2641                 } while(x < xe);
2642         {
2643                 if (b->wds >= b->maxwds) {
2644                         b1 = Balloc(b->k+1 MTa);
2645                         Bcopy(b1,b);
2646                         Bfree(b MTa);
2647                         b = b1;
2648                         }
2649                 b->x[b->wds++] = 1;
2650                 }
2651         return b;
2652         }
2653
2654 #endif /*}*/
2655
2656 #ifndef NO_HEX_FP /*{*/
2657
2658  static void
2659 rshift(Bigint *b, int k)
2660 {
2661         ULong *x, *x1, *xe, y;
2662         int n;
2663
2664         x = x1 = b->x;
2665         n = k >> kshift;
2666         if (n < b->wds) {
2667                 xe = x + b->wds;
2668                 x += n;
2669                 if (k &= kmask) {
2670                         n = 32 - k;
2671                         y = *x++ >> k;
2672                         while(x < xe) {
2673                                 *x1++ = (y | (*x << n)) & 0xffffffff;
2674                                 y = *x++ >> k;
2675                                 }
2676                         if ((*x1 = y) !=0)
2677                                 x1++;
2678                         }
2679                 else
2680                         while(x < xe)
2681                                 *x1++ = *x++;
2682                 }
2683         if ((b->wds = x1 - b->x) == 0)
2684                 b->x[0] = 0;
2685         }
2686
2687  static ULong
2688 any_on(Bigint *b, int k)
2689 {
2690         int n, nwds;
2691         ULong *x, *x0, x1, x2;
2692
2693         x = b->x;
2694         nwds = b->wds;
2695         n = k >> kshift;
2696         if (n > nwds)
2697                 n = nwds;
2698         else if (n < nwds && (k &= kmask)) {
2699                 x1 = x2 = x[n];
2700                 x1 >>= k;
2701                 x1 <<= k;
2702                 if (x1 != x2)
2703                         return 1;
2704                 }
2705         x0 = x;
2706         x += n;
2707         while(x > x0)
2708                 if (*--x)
2709                         return 1;
2710         return 0;
2711         }
2712
2713 enum {  /* rounding values: same as FLT_ROUNDS */
2714         Round_zero = 0,
2715         Round_near = 1,
2716         Round_up = 2,
2717         Round_down = 3
2718         };
2719
2720  void
2721 gethex( const char **sp, U *rvp, int rounding, int sign MTd)
2722 {
2723         Bigint *b;
2724         const unsigned char *decpt, *s0, *s, *s1;
2725         Long e, e1;
2726         ULong L, lostbits, *x;
2727         int big, denorm, esign, havedig, k, n, nbits, up, zret;
2728 #ifdef IBM
2729         int j;
2730 #endif
2731         enum {
2732 #ifdef IEEE_Arith /*{{*/
2733                 emax = 0x7fe - Bias - P + 1,
2734                 emin = Emin - P + 1
2735 #else /*}{*/
2736                 emin = Emin - P,
2737 #ifdef VAX
2738                 emax = 0x7ff - Bias - P + 1
2739 #endif
2740 #ifdef IBM
2741                 emax = 0x7f - Bias - P
2742 #endif
2743 #endif /*}}*/
2744                 };
2745 #ifdef USE_LOCALE
2746         int i;
2747 #ifdef NO_LOCALE_CACHE
2748         const unsigned char *decimalpoint = (unsigned char*)
2749                 localeconv()->decimal_point;
2750 #else
2751         const unsigned char *decimalpoint;
2752         static unsigned char *decimalpoint_cache;
2753         if (!(s0 = decimalpoint_cache)) {
2754                 s0 = (unsigned char*)localeconv()->decimal_point;
2755                 if ((decimalpoint_cache = (unsigned char*)
2756                                 MALLOC(strlen((const char*)s0) + 1))) {
2757                         strcpy((char*)decimalpoint_cache, (const char*)s0);
2758                         s0 = decimalpoint_cache;
2759                         }
2760                 }
2761         decimalpoint = s0;
2762 #endif
2763 #endif
2764
2765         /**** if (!hexdig['0']) hexdig_init(); ****/
2766         havedig = 0;
2767         s0 = *(const unsigned char **)sp + 2;
2768         while(s0[havedig] == '0')
2769                 havedig++;
2770         s0 += havedig;
2771         s = s0;
2772         decpt = 0;
2773         zret = 0;
2774         e = 0;
2775         if (hexdig[*s])
2776                 havedig++;
2777         else {
2778                 zret = 1;
2779 #ifdef USE_LOCALE
2780                 for(i = 0; decimalpoint[i]; ++i) {
2781                         if (s[i] != decimalpoint[i])
2782                                 goto pcheck;
2783                         }
2784                 decpt = s += i;
2785 #else
2786                 if (*s != '.')
2787                         goto pcheck;
2788                 decpt = ++s;
2789 #endif
2790                 if (!hexdig[*s])
2791                         goto pcheck;
2792                 while(*s == '0')
2793                         s++;
2794                 if (hexdig[*s])
2795                         zret = 0;
2796                 havedig = 1;
2797                 s0 = s;
2798                 }
2799         while(hexdig[*s])
2800                 s++;
2801 #ifdef USE_LOCALE
2802         if (*s == *decimalpoint && !decpt) {
2803                 for(i = 1; decimalpoint[i]; ++i) {
2804                         if (s[i] != decimalpoint[i])
2805                                 goto pcheck;
2806                         }
2807                 decpt = s += i;
2808 #else
2809         if (*s == '.' && !decpt) {
2810                 decpt = ++s;
2811 #endif
2812                 while(hexdig[*s])
2813                         s++;
2814                 }/*}*/
2815         if (decpt)
2816                 e = -(((Long)(s-decpt)) << 2);
2817  pcheck:
2818         s1 = s;
2819         big = esign = 0;
2820         switch(*s) {
2821           case 'p':
2822           case 'P':
2823                 switch(*++s) {
2824                   case '-':
2825                         esign = 1;
2826                         /* no break */
2827                   case '+':
2828                         s++;
2829                   }
2830                 if ((n = hexdig[*s]) == 0 || n > 0x19) {
2831                         s = s1;
2832                         break;
2833                         }
2834                 e1 = n - 0x10;
2835                 while((n = hexdig[*++s]) !=0 && n <= 0x19) {
2836                         if (e1 & 0xf8000000)
2837                                 big = 1;
2838                         e1 = 10*e1 + n - 0x10;
2839                         }
2840                 if (esign)
2841                         e1 = -e1;
2842                 e += e1;
2843           }
2844         *sp = (char*)s;
2845         if (!havedig)
2846                 *sp = (char*)s0 - 1;
2847         if (zret)
2848                 goto retz1;
2849         if (big) {
2850                 if (esign) {
2851 #ifdef IEEE_Arith
2852                         switch(rounding) {
2853                           case Round_up:
2854                                 if (sign)
2855                                         break;
2856                                 goto ret_tiny;
2857                           case Round_down:
2858                                 if (!sign)
2859                                         break;
2860                                 goto ret_tiny;
2861                           }
2862 #endif
2863                         goto retz;
2864 #ifdef IEEE_Arith
2865  ret_tinyf:
2866                         Bfree(b MTa);
2867  ret_tiny:
2868                         Set_errno(ERANGE);
2869                         word0(rvp) = 0;
2870                         word1(rvp) = 1;
2871                         return;
2872 #endif /* IEEE_Arith */
2873                         }
2874                 switch(rounding) {
2875                   case Round_near:
2876                         goto ovfl1;
2877                   case Round_up:
2878                         if (!sign)
2879                                 goto ovfl1;
2880                         goto ret_big;
2881                   case Round_down:
2882                         if (sign)
2883                                 goto ovfl1;
2884                         goto ret_big;
2885                   }
2886  ret_big:
2887                 word0(rvp) = Big0;
2888                 word1(rvp) = Big1;
2889                 return;
2890                 }
2891         n = s1 - s0 - 1;
2892         for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
2893                 k++;
2894         b = Balloc(k MTa);
2895         x = b->x;
2896         n = 0;
2897         L = 0;
2898 #ifdef USE_LOCALE
2899         for(i = 0; decimalpoint[i+1]; ++i);
2900 #endif
2901         while(s1 > s0) {
2902 #ifdef USE_LOCALE
2903                 if (*--s1 == decimalpoint[i]) {
2904                         s1 -= i;
2905                         continue;
2906                         }
2907 #else
2908                 if (*--s1 == '.')
2909                         continue;
2910 #endif
2911                 if (n == ULbits) {
2912                         *x++ = L;
2913                         L = 0;
2914                         n = 0;
2915                         }
2916                 L |= (hexdig[*s1] & 0x0f) << n;
2917                 n += 4;
2918                 }
2919         *x++ = L;
2920         b->wds = n = x - b->x;
2921         n = ULbits*n - hi0bits(L);
2922         nbits = Nbits;
2923         lostbits = 0;
2924         x = b->x;
2925         if (n > nbits) {
2926                 n -= nbits;
2927                 if (any_on(b,n)) {
2928                         lostbits = 1;
2929                         k = n - 1;
2930                         if (x[k>>kshift] & 1 << (k & kmask)) {
2931                                 lostbits = 2;
2932                                 if (k > 0 && any_on(b,k))
2933                                         lostbits = 3;
2934                                 }
2935                         }
2936                 rshift(b, n);
2937                 e += n;
2938                 }
2939         else if (n < nbits) {
2940                 n = nbits - n;
2941                 b = lshift(b, n MTa);
2942                 e -= n;
2943                 x = b->x;
2944                 }
2945         if (e > emax) {
2946  ovfl:
2947                 Bfree(b MTa);
2948  ovfl1:
2949                 Set_errno(ERANGE);
2950 #ifdef Honor_FLT_ROUNDS
2951                 switch (rounding) {
2952                   case Round_zero:
2953                         goto ret_big;
2954                   case Round_down:
2955                         if (!sign)
2956                                 goto ret_big;
2957                         break;
2958                   case Round_up:
2959                         if (sign)
2960                                 goto ret_big;
2961                   }
2962 #endif
2963                 word0(rvp) = Exp_mask;
2964                 word1(rvp) = 0;
2965                 return;
2966                 }
2967         denorm = 0;
2968         if (e < emin) {
2969                 denorm = 1;
2970                 n = emin - e;
2971                 if (n >= nbits) {
2972 #ifdef IEEE_Arith /*{*/
2973                         switch (rounding) {
2974                           case Round_near:
2975                                 if (n == nbits && (n < 2 || lostbits || any_on(b,n-1)))
2976                                         goto ret_tinyf;
2977                                 break;
2978                           case Round_up:
2979                                 if (!sign)
2980                                         goto ret_tinyf;
2981                                 break;
2982                           case Round_down:
2983                                 if (sign)
2984                                         goto ret_tinyf;
2985                           }
2986 #endif /* } IEEE_Arith */
2987                         Bfree(b MTa);
2988  retz:
2989                         Set_errno(ERANGE);
2990  retz1:
2991                         rvp->d = 0.;
2992                         return;
2993                         }
2994                 k = n - 1;
2995                 if (lostbits)
2996                         lostbits = 1;
2997                 else if (k > 0)
2998                         lostbits = any_on(b,k);
2999                 if (x[k>>kshift] & 1 << (k & kmask))
3000                         lostbits |= 2;
3001                 nbits -= n;
3002                 rshift(b,n);
3003                 e = emin;
3004                 }
3005         if (lostbits) {
3006                 up = 0;
3007                 switch(rounding) {
3008                   case Round_zero:
3009                         break;
3010                   case Round_near:
3011                         if (lostbits & 2
3012                          && (lostbits & 1) | (x[0] & 1))
3013                                 up = 1;
3014                         break;
3015                   case Round_up:
3016                         up = 1 - sign;
3017                         break;
3018                   case Round_down:
3019                         up = sign;
3020                   }
3021                 if (up) {
3022                         k = b->wds;
3023                         b = increment(b MTa);
3024                         x = b->x;
3025                         if (denorm) {
3026 #if 0
3027                                 if (nbits == Nbits - 1
3028                                  && x[nbits >> kshift] & 1 << (nbits & kmask))
3029                                         denorm = 0; /* not currently used */
3030 #endif
3031                                 }
3032                         else if (b->wds > k
3033                          || ((n = nbits & kmask) !=0
3034                              && hi0bits(x[k-1]) < 32-n)) {
3035                                 rshift(b,1);
3036                                 if (++e > Emax)
3037                                         goto ovfl;
3038                                 }
3039                         }
3040                 }
3041 #ifdef IEEE_Arith
3042         if (denorm)
3043                 word0(rvp) = b->wds > 1 ? b->x[1] & ~0x100000 : 0;
3044         else
3045                 word0(rvp) = (b->x[1] & ~0x100000) | ((e + 0x3ff + 52) << 20);
3046         word1(rvp) = b->x[0];
3047 #endif
3048 #ifdef IBM
3049         if ((j = e & 3)) {
3050                 k = b->x[0] & ((1 << j) - 1);
3051                 rshift(b,j);
3052                 if (k) {
3053                         switch(rounding) {
3054                           case Round_up:
3055                                 if (!sign)
3056                                         increment(b);
3057                                 break;
3058                           case Round_down:
3059                                 if (sign)
3060                                         increment(b);
3061                                 break;
3062                           case Round_near:
3063                                 j = 1 << (j-1);
3064                                 if (k & j && ((k & (j-1)) | lostbits))
3065                                         increment(b);
3066                           }
3067                         }
3068                 }
3069         e >>= 2;
3070         word0(rvp) = b->x[1] | ((e + 65 + 13) << 24);
3071         word1(rvp) = b->x[0];
3072 #endif
3073 #ifdef VAX
3074         /* The next two lines ignore swap of low- and high-order 2 bytes. */
3075         /* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */
3076         /* word1(rvp) = b->x[0]; */
3077         word0(rvp) = ((b->x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b->x[1] << 16);
3078         word1(rvp) = (b->x[0] >> 16) | (b->x[0] << 16);
3079 #endif
3080         Bfree(b MTa);
3081         }
3082 #endif /*!NO_HEX_FP}*/
3083
3084  static int
3085 dshift(Bigint *b, int p2)
3086 {
3087         int rv = hi0bits(b->x[b->wds-1]) - 4;
3088         if (p2 > 0)
3089                 rv -= p2;
3090         return rv & kmask;
3091         }
3092
3093  static int
3094 quorem(Bigint *b, Bigint *S)
3095 {
3096         int n;
3097         ULong *bx, *bxe, q, *sx, *sxe;
3098 #ifdef ULLong
3099         ULLong borrow, carry, y, ys;
3100 #else
3101         ULong borrow, carry, y, ys;
3102 #ifdef Pack_32
3103         ULong si, z, zs;
3104 #endif
3105 #endif
3106
3107         n = S->wds;
3108 #ifdef DEBUG
3109         /*debug*/ if (b->wds > n)
3110         /*debug*/       Bug("oversize b in quorem");
3111 #endif
3112         if (b->wds < n)
3113                 return 0;
3114         sx = S->x;
3115         sxe = sx + --n;
3116         bx = b->x;
3117         bxe = bx + n;
3118         q = *bxe / (*sxe + 1);  /* ensure q <= true quotient */
3119 #ifdef DEBUG
3120 #ifdef NO_STRTOD_BIGCOMP
3121         /*debug*/ if (q > 9)
3122 #else
3123         /* An oversized q is possible when quorem is called from bigcomp and */
3124         /* the input is near, e.g., twice the smallest denormalized number. */
3125         /*debug*/ if (q > 15)
3126 #endif
3127         /*debug*/       Bug("oversized quotient in quorem");
3128 #endif
3129         if (q) {
3130                 borrow = 0;
3131                 carry = 0;
3132                 do {
3133 #ifdef ULLong
3134                         ys = *sx++ * (ULLong)q + carry;
3135                         carry = ys >> 32;
3136                         y = *bx - (ys & FFFFFFFF) - borrow;
3137                         borrow = y >> 32 & (ULong)1;
3138                         *bx++ = y & FFFFFFFF;
3139 #else
3140 #ifdef Pack_32
3141                         si = *sx++;
3142                         ys = (si & 0xffff) * q + carry;
3143                         zs = (si >> 16) * q + (ys >> 16);
3144                         carry = zs >> 16;
3145                         y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
3146                         borrow = (y & 0x10000) >> 16;
3147                         z = (*bx >> 16) - (zs & 0xffff) - borrow;
3148                         borrow = (z & 0x10000) >> 16;
3149                         Storeinc(bx, z, y);
3150 #else
3151                         ys = *sx++ * q + carry;
3152                         carry = ys >> 16;
3153                         y = *bx - (ys & 0xffff) - borrow;
3154                         borrow = (y & 0x10000) >> 16;
3155                         *bx++ = y & 0xffff;
3156 #endif
3157 #endif
3158                         }
3159                         while(sx <= sxe);
3160                 if (!*bxe) {
3161                         bx = b->x;
3162                         while(--bxe > bx && !*bxe)
3163                                 --n;
3164                         b->wds = n;
3165                         }
3166                 }
3167         if (cmp(b, S) >= 0) {
3168                 q++;
3169                 borrow = 0;
3170                 carry = 0;
3171                 bx = b->x;
3172                 sx = S->x;
3173                 do {
3174 #ifdef ULLong
3175                         ys = *sx++ + carry;
3176                         carry = ys >> 32;
3177                         y = *bx - (ys & FFFFFFFF) - borrow;
3178                         borrow = y >> 32 & (ULong)1;
3179                         *bx++ = y & FFFFFFFF;
3180 #else
3181 #ifdef Pack_32
3182                         si = *sx++;
3183                         ys = (si & 0xffff) + carry;
3184                         zs = (si >> 16) + (ys >> 16);
3185                         carry = zs >> 16;
3186                         y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
3187                         borrow = (y & 0x10000) >> 16;
3188                         z = (*bx >> 16) - (zs & 0xffff) - borrow;
3189                         borrow = (z & 0x10000) >> 16;
3190                         Storeinc(bx, z, y);
3191 #else
3192                         ys = *sx++ + carry;
3193                         carry = ys >> 16;
3194                         y = *bx - (ys & 0xffff) - borrow;
3195                         borrow = (y & 0x10000) >> 16;
3196                         *bx++ = y & 0xffff;
3197 #endif
3198 #endif
3199                         }
3200                         while(sx <= sxe);
3201                 bx = b->x;
3202                 bxe = bx + n;
3203                 if (!*bxe) {
3204                         while(--bxe > bx && !*bxe)
3205                                 --n;
3206                         b->wds = n;
3207                         }
3208                 }
3209         return q;
3210         }
3211
3212 #if defined(Avoid_Underflow) || !defined(NO_STRTOD_BIGCOMP) /*{*/
3213  static double
3214 sulp(U *x, BCinfo *bc)
3215 {
3216         U u;
3217         double rv;
3218         int i;
3219
3220         rv = ulp(x);
3221         if (!bc->scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0)
3222                 return rv; /* Is there an example where i <= 0 ? */
3223         word0(&u) = Exp_1 + (i << Exp_shift);
3224         word1(&u) = 0;
3225         return rv * u.d;
3226         }
3227 #endif /*}*/
3228
3229 #ifndef NO_STRTOD_BIGCOMP
3230  static void
3231 bigcomp(U *rv, const char *s0, BCinfo *bc MTd)
3232 {
3233         Bigint *b, *d;
3234         int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase;
3235
3236         dsign = bc->dsign;
3237         nd = bc->nd;
3238         nd0 = bc->nd0;
3239         p5 = nd + bc->e0 - 1;
3240         speccase = 0;
3241 #ifndef Sudden_Underflow
3242         if (rv->d == 0.) {      /* special case: value near underflow-to-zero */
3243                                 /* threshold was rounded to zero */
3244                 b = i2b(1 MTa);
3245                 p2 = Emin - P + 1;
3246                 bbits = 1;
3247 #ifdef Avoid_Underflow
3248                 word0(rv) = (P+2) << Exp_shift;
3249 #else
3250                 word1(rv) = 1;
3251 #endif
3252                 i = 0;
3253 #ifdef Honor_FLT_ROUNDS
3254                 if (bc->rounding == 1)
3255 #endif
3256                         {
3257                         speccase = 1;
3258                         --p2;
3259                         dsign = 0;
3260                         goto have_i;
3261                         }
3262                 }
3263         else
3264 #endif
3265                 b = d2b(rv, &p2, &bbits MTa);
3266 #ifdef Avoid_Underflow
3267         p2 -= bc->scale;
3268 #endif
3269         /* floor(log2(rv)) == bbits - 1 + p2 */
3270         /* Check for denormal case. */
3271         i = P - bbits;
3272         if (i > (j = P - Emin - 1 + p2)) {
3273 #ifdef Sudden_Underflow
3274                 Bfree(b MTa);
3275                 b = i2b(1 MTa);
3276                 p2 = Emin;
3277                 i = P - 1;
3278 #ifdef Avoid_Underflow
3279                 word0(rv) = (1 + bc->scale) << Exp_shift;
3280 #else
3281                 word0(rv) = Exp_msk1;
3282 #endif
3283                 word1(rv) = 0;
3284 #else
3285                 i = j;
3286 #endif
3287                 }
3288 #ifdef Honor_FLT_ROUNDS
3289         if (bc->rounding != 1) {
3290                 if (i > 0)
3291                         b = lshift(b, i MTa);
3292                 if (dsign)
3293                         b = increment(b MTa);
3294                 }
3295         else
3296 #endif
3297                 {
3298                 b = lshift(b, ++i MTa);
3299                 b->x[0] |= 1;
3300                 }
3301 #ifndef Sudden_Underflow
3302  have_i:
3303 #endif
3304         p2 -= p5 + i;
3305         d = i2b(1 MTa);
3306         /* Arrange for convenient computation of quotients:
3307          * shift left if necessary so divisor has 4 leading 0 bits.
3308          */
3309         if (p5 > 0)
3310                 d = pow5mult(d, p5 MTa);
3311         else if (p5 < 0)
3312                 b = pow5mult(b, -p5 MTa);
3313         if (p2 > 0) {
3314                 b2 = p2;
3315                 d2 = 0;
3316                 }
3317         else {
3318                 b2 = 0;
3319                 d2 = -p2;
3320                 }
3321         i = dshift(d, d2);
3322         if ((b2 += i) > 0)
3323                 b = lshift(b, b2 MTa);
3324         if ((d2 += i) > 0)
3325                 d = lshift(d, d2 MTa);
3326
3327         /* Now b/d = exactly half-way between the two floating-point values */
3328         /* on either side of the input string.  Compute first digit of b/d. */
3329
3330         if (!(dig = quorem(b,d))) {
3331                 b = multadd(b, 10, 0 MTa);      /* very unlikely */
3332                 dig = quorem(b,d);
3333                 }
3334
3335         /* Compare b/d with s0 */
3336
3337         for(i = 0; i < nd0; ) {
3338                 if ((dd = s0[i++] - '0' - dig))
3339                         goto ret;
3340                 if (!b->x[0] && b->wds == 1) {
3341                         if (i < nd)
3342                                 dd = 1;
3343                         goto ret;
3344                         }
3345                 b = multadd(b, 10, 0 MTa);
3346                 dig = quorem(b,d);
3347                 }
3348         for(j = bc->dp1; i++ < nd;) {
3349                 if ((dd = s0[j++] - '0' - dig))
3350                         goto ret;
3351                 if (!b->x[0] && b->wds == 1) {
3352                         if (i < nd)
3353                                 dd = 1;
3354                         goto ret;
3355                         }
3356                 b = multadd(b, 10, 0 MTa);
3357                 dig = quorem(b,d);
3358                 }
3359         if (dig > 0 || b->x[0] || b->wds > 1)
3360                 dd = -1;
3361  ret:
3362         Bfree(b MTa);
3363         Bfree(d MTa);
3364 #ifdef Honor_FLT_ROUNDS
3365         if (bc->rounding != 1) {
3366                 if (dd < 0) {
3367                         if (bc->rounding == 0) {
3368                                 if (!dsign)
3369                                         goto retlow1;
3370                                 }
3371                         else if (dsign)
3372                                 goto rethi1;
3373                         }
3374                 else if (dd > 0) {
3375                         if (bc->rounding == 0) {
3376                                 if (dsign)
3377                                         goto rethi1;
3378                                 goto ret1;
3379                                 }
3380                         if (!dsign)
3381                                 goto rethi1;
3382                         dval(rv) += 2.*sulp(rv,bc);
3383                         }
3384                 else {
3385                         bc->inexact = 0;
3386                         if (dsign)
3387                                 goto rethi1;
3388                         }
3389                 }
3390         else
3391 #endif
3392         if (speccase) {
3393                 if (dd <= 0)
3394                         rv->d = 0.;
3395                 }
3396         else if (dd < 0) {
3397                 if (!dsign)     /* does not happen for round-near */
3398 retlow1:
3399                         dval(rv) -= sulp(rv,bc);
3400                 }
3401         else if (dd > 0) {
3402                 if (dsign) {
3403  rethi1:
3404                         dval(rv) += sulp(rv,bc);
3405                         }
3406                 }
3407         else {
3408                 /* Exact half-way case:  apply round-even rule. */
3409                 if ((j = ((word0(rv) & Exp_mask) >> Exp_shift) - bc->scale) <= 0) {
3410                         i = 1 - j;
3411                         if (i <= 31) {
3412                                 if (word1(rv) & (0x1 << i))
3413                                         goto odd;
3414                                 }
3415                         else if (word0(rv) & (0x1 << (i-32)))
3416                                 goto odd;
3417                         }
3418                 else if (word1(rv) & 1) {
3419  odd:
3420                         if (dsign)
3421                                 goto rethi1;
3422                         goto retlow1;
3423                         }
3424                 }
3425
3426 #ifdef Honor_FLT_ROUNDS
3427  ret1:
3428 #endif
3429         return;
3430         }
3431 #endif /* NO_STRTOD_BIGCOMP */
3432
3433  double
3434 strtod(const char *s00, char **se)
3435 {
3436         int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1;
3437         int esign, i, j, k, nd, nd0, nf, nz, nz0, nz1, sign;
3438         const char *s, *s0, *s1;
3439         double aadj, aadj1;
3440         Long L;
3441         U aadj2, adj, rv, rv0;
3442         ULong y, z;
3443         BCinfo bc;
3444         Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
3445 #ifdef USE_BF96
3446         ULLong bhi, blo, brv, t00, t01, t02, t10, t11, terv, tg, tlo, yz;
3447         const BF96 *p10;
3448         int bexact, erv;
3449 #endif
3450 #ifdef Avoid_Underflow
3451         ULong Lsb, Lsb1;
3452 #endif
3453 #ifdef SET_INEXACT
3454         int oldinexact;
3455 #endif
3456 #ifndef NO_STRTOD_BIGCOMP
3457         int req_bigcomp = 0;
3458 #endif
3459 #ifdef MULTIPLE_THREADS
3460         ThInfo *TI = 0;
3461 #endif
3462 #ifdef Honor_FLT_ROUNDS /*{*/
3463 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
3464         bc.rounding = Flt_Rounds;
3465 #else /*}{*/
3466         bc.rounding = 1;
3467         switch(fegetround()) {
3468           case FE_TOWARDZERO:   bc.rounding = 0; break;
3469           case FE_UPWARD:       bc.rounding = 2; break;
3470           case FE_DOWNWARD:     bc.rounding = 3;
3471           }
3472 #endif /*}}*/
3473 #endif /*}*/
3474 #ifdef USE_LOCALE
3475         const char *s2;
3476 #endif
3477
3478         sign = nz0 = nz1 = nz = bc.dplen = bc.uflchk = 0;
3479         dval(&rv) = 0.;
3480         for(s = s00;;s++) switch(*s) {
3481                 case '-':
3482                         sign = 1;
3483                         /* no break */
3484                 case '+':
3485                         if (*++s)
3486                                 goto break2;
3487                         /* no break */
3488                 case 0:
3489                         goto ret0;
3490                 case '\t':
3491                 case '\n':
3492                 case '\v':
3493                 case '\f':
3494                 case '\r':
3495                 case ' ':
3496                         continue;
3497                 default:
3498                         goto break2;
3499                 }
3500  break2:
3501         if (*s == '0') {
3502 #ifndef NO_HEX_FP /*{*/
3503                 switch(s[1]) {
3504                   case 'x':
3505                   case 'X':
3506 #ifdef Honor_FLT_ROUNDS
3507                         gethex(&s, &rv, bc.rounding, sign MTb);
3508 #else
3509                         gethex(&s, &rv, 1, sign MTb);
3510 #endif
3511                         goto ret;
3512                   }
3513 #endif /*}*/
3514                 nz0 = 1;
3515                 while(*++s == '0') ;
3516                 if (!*s)
3517                         goto ret;
3518                 }
3519         s0 = s;
3520         nd = nf = 0;
3521 #ifdef USE_BF96
3522         yz = 0;
3523         for(; (c = *s) >= '0' && c <= '9'; nd++, s++)
3524                 if (nd < 19)
3525                         yz = 10*yz + c - '0';
3526 #else
3527         y = z = 0;
3528         for(; (c = *s) >= '0' && c <= '9'; nd++, s++)
3529                 if (nd < 9)
3530                         y = 10*y + c - '0';
3531                 else if (nd < DBL_DIG + 2)
3532                         z = 10*z + c - '0';
3533 #endif
3534         nd0 = nd;
3535         bc.dp0 = bc.dp1 = s - s0;
3536         for(s1 = s; s1 > s0 && *--s1 == '0'; )
3537                 ++nz1;
3538 #ifdef USE_LOCALE
3539         s1 = localeconv()->decimal_point;
3540         if (c == *s1) {
3541                 c = '.';
3542                 if (*++s1) {
3543                         s2 = s;
3544                         for(;;) {
3545                                 if (*++s2 != *s1) {
3546                                         c = 0;
3547                                         break;
3548                                         }
3549                                 if (!*++s1) {
3550                                         s = s2;
3551                                         break;
3552                                         }
3553                                 }
3554                         }
3555                 }
3556 #endif
3557         if (c == '.') {
3558                 c = *++s;
3559                 bc.dp1 = s - s0;
3560                 bc.dplen = bc.dp1 - bc.dp0;
3561                 if (!nd) {
3562                         for(; c == '0'; c = *++s)
3563                                 nz++;
3564                         if (c > '0' && c <= '9') {
3565                                 bc.dp0 = s0 - s;
3566                                 bc.dp1 = bc.dp0 + bc.dplen;
3567                                 s0 = s;
3568                                 nf += nz;
3569                                 nz = 0;
3570                                 goto have_dig;
3571                                 }
3572                         goto dig_done;
3573                         }
3574                 for(; c >= '0' && c <= '9'; c = *++s) {
3575  have_dig:
3576                         nz++;
3577                         if (c -= '0') {
3578                                 nf += nz;
3579                                 i = 1;
3580 #ifdef USE_BF96
3581                                 for(; i < nz; ++i) {
3582                                         if (++nd <= 19)
3583                                                 yz *= 10;
3584                                         }
3585                                 if (++nd <= 19)
3586                                         yz = 10*yz + c;
3587 #else
3588                                 for(; i < nz; ++i) {
3589                                         if (nd++ < 9)
3590                                                 y *= 10;
3591                                         else if (nd <= DBL_DIG + 2)
3592                                                 z *= 10;
3593                                         }
3594                                 if (nd++ < 9)
3595                                         y = 10*y + c;
3596                                 else if (nd <= DBL_DIG + 2)
3597                                         z = 10*z + c;
3598 #endif
3599                                 nz = nz1 = 0;
3600                                 }
3601                         }
3602                 }
3603  dig_done:
3604         e = 0;
3605         if (c == 'e' || c == 'E') {
3606                 if (!nd && !nz && !nz0) {
3607                         goto ret0;
3608                         }
3609                 s00 = s;
3610                 esign = 0;
3611                 switch(c = *++s) {
3612                         case '-':
3613                                 esign = 1;
3614                         case '+':
3615                                 c = *++s;
3616                         }
3617                 if (c >= '0' && c <= '9') {
3618                         while(c == '0')
3619                                 c = *++s;
3620                         if (c > '0' && c <= '9') {
3621                                 L = c - '0';
3622                                 s1 = s;
3623                                 while((c = *++s) >= '0' && c <= '9')
3624                                         L = 10*L + c - '0';
3625                                 if (s - s1 > 8 || L > 19999)
3626                                         /* Avoid confusion from exponents
3627                                          * so large that e might overflow.
3628                                          */
3629                                         e = 19999; /* safe for 16 bit ints */
3630                                 else
3631                                         e = (int)L;
3632                                 if (esign)
3633                                         e = -e;
3634                                 }
3635                         else
3636                                 e = 0;
3637                         }
3638                 else
3639                         s = s00;
3640                 }
3641         if (!nd) {
3642                 if (!nz && !nz0) {
3643 #ifdef INFNAN_CHECK /*{*/
3644                         /* Check for Nan and Infinity */
3645                         if (!bc.dplen)
3646                          switch(c) {
3647                           case 'i':
3648                           case 'I':
3649                                 if (match(&s,"nf")) {
3650                                         --s;
3651                                         if (!match(&s,"inity"))
3652                                                 ++s;
3653                                         word0(&rv) = 0x7ff00000;
3654                                         word1(&rv) = 0;
3655                                         goto ret;
3656                                         }
3657                                 break;
3658                           case 'n':
3659                           case 'N':
3660                                 if (match(&s, "an")) {
3661                                         word0(&rv) = NAN_WORD0;
3662                                         word1(&rv) = NAN_WORD1;
3663 #ifndef No_Hex_NaN
3664                                         if (*s == '(') /*)*/
3665                                                 hexnan(&rv, &s);
3666 #endif
3667                                         goto ret;
3668                                         }
3669                           }
3670 #endif /*} INFNAN_CHECK */
3671  ret0:
3672                         s = s00;
3673                         sign = 0;
3674                         }
3675                 goto ret;
3676                 }
3677         bc.e0 = e1 = e -= nf;
3678
3679         /* Now we have nd0 digits, starting at s0, followed by a
3680          * decimal point, followed by nd-nd0 digits.  The number we're
3681          * after is the integer represented by those digits times
3682          * 10**e */
3683
3684         if (!nd0)
3685                 nd0 = nd;
3686 #ifndef USE_BF96
3687         k = nd < DBL_DIG + 2 ? nd : DBL_DIG + 2;
3688         dval(&rv) = y;
3689         if (k > 9) {
3690 #ifdef SET_INEXACT
3691                 if (k > DBL_DIG)
3692                         oldinexact = get_inexact();
3693 #endif
3694                 dval(&rv) = tens[k - 9] * dval(&rv) + z;
3695                 }
3696 #endif
3697         bd0 = 0;
3698         if (nd <= DBL_DIG
3699 #ifndef RND_PRODQUOT
3700 #ifndef Honor_FLT_ROUNDS
3701                 && Flt_Rounds == 1
3702 #endif
3703 #endif
3704                         ) {
3705 #ifdef USE_BF96
3706                 dval(&rv) = yz;
3707 #endif
3708                 if (!e)
3709                         goto ret;
3710 #ifndef ROUND_BIASED_without_Round_Up
3711                 if (e > 0) {
3712                         if (e <= Ten_pmax) {
3713 #ifdef SET_INEXACT
3714                                 bc.inexact = 0;
3715                                 oldinexact = 1;
3716 #endif
3717 #ifdef VAX
3718                                 goto vax_ovfl_check;
3719 #else
3720 #ifdef Honor_FLT_ROUNDS
3721                                 /* round correctly FLT_ROUNDS = 2 or 3 */
3722                                 if (sign) {
3723                                         rv.d = -rv.d;
3724                                         sign = 0;
3725                                         }
3726 #endif
3727                                 /* rv = */ rounded_product(dval(&rv), tens[e]);
3728                                 goto ret;
3729 #endif
3730                                 }
3731                         i = DBL_DIG - nd;
3732                         if (e <= Ten_pmax + i) {
3733                                 /* A fancier test would sometimes let us do
3734                                  * this for larger i values.
3735                                  */
3736 #ifdef SET_INEXACT
3737                                 bc.inexact = 0;
3738                                 oldinexact = 1;
3739 #endif
3740 #ifdef Honor_FLT_ROUNDS
3741                                 /* round correctly FLT_ROUNDS = 2 or 3 */
3742                                 if (sign) {
3743                                         rv.d = -rv.d;
3744                                         sign = 0;
3745                                         }
3746 #endif
3747                                 e -= i;
3748                                 dval(&rv) *= tens[i];
3749 #ifdef VAX
3750                                 /* VAX exponent range is so narrow we must
3751                                  * worry about overflow here...
3752                                  */
3753  vax_ovfl_check:
3754                                 word0(&rv) -= P*Exp_msk1;
3755                                 /* rv = */ rounded_product(dval(&rv), tens[e]);
3756                                 if ((word0(&rv) & Exp_mask)
3757                                  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
3758                                         goto ovfl;
3759                                 word0(&rv) += P*Exp_msk1;
3760 #else
3761                                 /* rv = */ rounded_product(dval(&rv), tens[e]);
3762 #endif
3763                                 goto ret;
3764                                 }
3765                         }
3766 #ifndef Inaccurate_Divide
3767                 else if (e >= -Ten_pmax) {
3768 #ifdef SET_INEXACT
3769                                 bc.inexact = 0;
3770                                 oldinexact = 1;
3771 #endif
3772 #ifdef Honor_FLT_ROUNDS
3773                         /* round correctly FLT_ROUNDS = 2 or 3 */
3774                         if (sign) {
3775                                 rv.d = -rv.d;
3776                                 sign = 0;
3777                                 }
3778 #endif
3779                         /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
3780                         goto ret;
3781                         }
3782 #endif
3783 #endif /* ROUND_BIASED_without_Round_Up */
3784                 }
3785 #ifdef USE_BF96
3786         k = nd < 19 ? nd : 19;
3787 #endif
3788         e1 += nd - k;   /* scale factor = 10^e1 */
3789
3790 #ifdef IEEE_Arith
3791 #ifdef SET_INEXACT
3792         bc.inexact = 1;
3793 #ifndef USE_BF96
3794         if (k <= DBL_DIG)
3795 #endif
3796                 oldinexact = get_inexact();
3797 #endif
3798 #ifdef Honor_FLT_ROUNDS
3799         if (bc.rounding >= 2) {
3800                 if (sign)
3801                         bc.rounding = bc.rounding == 2 ? 0 : 2;
3802                 else
3803                         if (bc.rounding != 2)
3804                                 bc.rounding = 0;
3805                 }
3806 #endif
3807 #endif /*IEEE_Arith*/
3808
3809 #ifdef USE_BF96 /*{*/
3810         Debug(++dtoa_stats[0]);
3811         i = e1 + 342;
3812         if (i < 0)
3813                 goto undfl;
3814         if (i > 650)
3815                 goto ovfl;
3816         p10 = &pten[i];
3817         brv = yz;
3818         /* shift brv left, with i =  number of bits shifted */
3819         i = 0;
3820         if (!(brv & 0xffffffff00000000ull)) {
3821                 i = 32;
3822                 brv <<= 32;
3823                 }
3824         if (!(brv & 0xffff000000000000ull)) {
3825                 i += 16;
3826                 brv <<= 16;
3827                 }
3828         if (!(brv & 0xff00000000000000ull)) {
3829                 i += 8;
3830                 brv <<= 8;
3831                 }
3832         if (!(brv & 0xf000000000000000ull)) {
3833                 i += 4;
3834                 brv <<= 4;
3835                 }
3836         if (!(brv & 0xc000000000000000ull)) {
3837                 i += 2;
3838                 brv <<= 2;
3839                 }
3840         if (!(brv & 0x8000000000000000ull)) {
3841                 i += 1;
3842                 brv <<= 1;
3843                 }
3844         erv = (64 + 0x3fe) + p10->e - i;
3845         if (erv <= 0 && nd > 19)
3846                 goto many_digits; /* denormal: may need to look at all digits */
3847         bhi = brv >> 32;
3848         blo = brv & 0xffffffffull;
3849         /* Unsigned 32-bit ints lie in [0,2^32-1] and */
3850         /* unsigned 64-bit ints lie in [0, 2^64-1].  The product of two unsigned */
3851         /* 32-bit ints is <= 2^64 - 2*2^32-1 + 1 = 2^64 - 1 - 2*(2^32 - 1), so */
3852         /* we can add two unsigned 32-bit ints to the product of two such ints, */
3853         /* and 64 bits suffice to contain the result. */
3854         t01 = bhi * p10->b1;
3855         t10 = blo * p10->b0 + (t01 & 0xffffffffull);
3856         t00 = bhi * p10->b0 + (t01 >> 32) + (t10 >> 32);
3857         if (t00 & 0x8000000000000000ull) {
3858                 if ((t00 & 0x3ff) && (~t00 & 0x3fe)) { /* unambiguous result? */
3859                         if (nd > 19 && ((t00 + (1<<i) + 2) & 0x400) ^ (t00 & 0x400))
3860                                 goto many_digits;
3861                         if (erv <= 0)
3862                                 goto denormal;
3863 #ifdef Honor_FLT_ROUNDS
3864                         switch(bc.rounding) {
3865                           case 0: goto noround;
3866                           case 2: goto roundup;
3867                           }
3868 #endif
3869                         if (t00 & 0x400 && t00 & 0xbff)
3870                                 goto roundup;
3871                         goto noround;
3872                         }
3873                 }
3874         else {
3875                 if ((t00 & 0x1ff) && (~t00 & 0x1fe)) { /* unambiguous result? */
3876                         if (nd > 19 && ((t00 + (1<<i) + 2) & 0x200) ^ (t00 & 0x200))
3877                                 goto many_digits;
3878                         if (erv <= 1)
3879                                 goto denormal1;
3880 #ifdef Honor_FLT_ROUNDS
3881                         switch(bc.rounding) {
3882                           case 0: goto noround1;
3883                           case 2: goto roundup1;
3884                           }
3885 #endif
3886                         if (t00 & 0x200)
3887                                 goto roundup1;
3888                         goto noround1;
3889                         }
3890                 }
3891         /* 3 multiplies did not suffice; try a 96-bit approximation */
3892         Debug(++dtoa_stats[1]);
3893         t02 = bhi * p10->b2;
3894         t11 = blo * p10->b1 + (t02 & 0xffffffffull);
3895         bexact = 1;
3896         if (e1 < 0 || e1 > 41 || (t10 | t11) & 0xffffffffull || nd > 19)
3897                 bexact = 0;
3898         tlo = (t10 & 0xffffffffull) + (t02 >> 32) + (t11 >> 32);
3899         if (!bexact && (tlo + 0x10) >> 32 > tlo >> 32)
3900                 goto many_digits;
3901         t00 += tlo >> 32;
3902         if (t00 & 0x8000000000000000ull) {
3903                 if (erv <= 0) { /* denormal result */
3904                         if (nd >= 20 || !((tlo & 0xfffffff0) | (t00 & 0x3ff)))
3905                                 goto many_digits;
3906  denormal:
3907                         if (erv <= -52) {
3908 #ifdef Honor_FLT_ROUNDS
3909                                 switch(bc.rounding) {
3910                                   case 0: goto undfl;
3911                                   case 2: goto tiniest;
3912                                   }
3913 #endif
3914                                 if (erv < -52 || !(t00 & 0x7fffffffffffffffull))
3915                                         goto undfl;
3916                                 goto tiniest;
3917                                 }
3918                         tg = 1ull << (11 - erv);
3919                         t00 &= ~(tg - 1); /* clear low bits */
3920 #ifdef Honor_FLT_ROUNDS
3921                         switch(bc.rounding) {
3922                           case 0: goto noround_den;
3923                           case 2: goto roundup_den;
3924                           }
3925 #endif
3926                         if (t00 & tg) {
3927 #ifdef Honor_FLT_ROUNDS
3928  roundup_den:
3929 #endif
3930                                 t00 += tg << 1;
3931                                 if (!(t00 & 0x8000000000000000ull)) {
3932                                         if (++erv > 0)
3933                                                 goto smallest_normal;
3934                                         t00 = 0x8000000000000000ull;
3935                                         }
3936                                 }
3937 #ifdef Honor_FLT_ROUNDS
3938  noround_den:
3939 #endif
3940                         LLval(&rv) = t00 >> (12 - erv);
3941                         Set_errno(ERANGE);
3942                         goto ret;
3943                         }
3944                 if (bexact) {
3945 #ifdef SET_INEXACT
3946                         if (!(t00 & 0x7ff) && !(tlo & 0xffffffffull)) {
3947                                 bc.inexact = 0;
3948                                 goto noround;
3949                                 }
3950 #endif
3951 #ifdef Honor_FLT_ROUNDS
3952                         switch(bc.rounding) {
3953                           case 2:
3954                                 if (t00 & 0x7ff)
3955                                         goto roundup;
3956                           case 0: goto noround;
3957                           }
3958 #endif
3959                         if (t00 & 0x400 && (tlo & 0xffffffff) | (t00 & 0xbff))
3960                                 goto roundup;
3961                         goto noround;
3962                         }
3963                 if ((tlo & 0xfffffff0) | (t00 & 0x3ff)
3964                  && (nd <= 19 ||  ((t00 + (1ull << i)) & 0xfffffffffffffc00ull)
3965                                 == (t00 & 0xfffffffffffffc00ull))) {
3966                         /* Unambiguous result. */
3967                         /* If nd > 19, then incrementing the 19th digit */
3968                         /* does not affect rv. */
3969 #ifdef Honor_FLT_ROUNDS
3970                         switch(bc.rounding) {
3971                           case 0: goto noround;
3972                           case 2: goto roundup;
3973                           }
3974 #endif
3975                         if (t00 & 0x400) { /* round up */
3976  roundup:
3977                                 t00 += 0x800;
3978                                 if (!(t00 & 0x8000000000000000ull)) {
3979                                         /* rounded up to a power of 2 */
3980                                         if (erv >= 0x7fe)
3981                                                 goto ovfl;
3982                                         terv = erv + 1;
3983                                         LLval(&rv) = terv << 52;
3984                                         goto ret;
3985                                         }
3986                                 }
3987  noround:
3988                         if (erv >= 0x7ff)
3989                                 goto ovfl;
3990                         terv = erv;
3991                         LLval(&rv) = (terv << 52) | ((t00 & 0x7ffffffffffff800ull) >> 11);
3992                         goto ret;
3993                         }
3994                 }
3995         else {
3996                 if (erv <= 1) { /* denormal result */
3997                         if (nd >= 20 || !((tlo & 0xfffffff0) | (t00 & 0x1ff)))
3998                                 goto many_digits;
3999  denormal1:
4000                         if (erv <= -51) {
4001 #ifdef Honor_FLT_ROUNDS
4002                                 switch(bc.rounding) {
4003                                   case 0: goto undfl;
4004                                   case 2: goto tiniest;
4005                                   }
4006 #endif
4007                                 if (erv < -51 || !(t00 & 0x3fffffffffffffffull))
4008                                         goto undfl;
4009  tiniest:
4010                                 LLval(&rv) = 1;
4011                                 Set_errno(ERANGE);
4012                                 goto ret;
4013                                 }
4014                         tg = 1ull << (11 - erv);
4015 #ifdef Honor_FLT_ROUNDS
4016                         switch(bc.rounding) {
4017                           case 0: goto noround1_den;
4018                           case 2: goto roundup1_den;
4019                           }
4020 #endif
4021                         if (t00 & tg) {
4022 #ifdef Honor_FLT_ROUNDS
4023  roundup1_den:
4024 #endif
4025                                 if (0x8000000000000000ull & (t00 += (tg<<1)) && erv == 1) {
4026
4027  smallest_normal:
4028                                         LLval(&rv) = 0x0010000000000000ull;
4029                                         goto ret;
4030                                         }
4031                                 }
4032 #ifdef Honor_FLT_ROUNDS
4033  noround1_den:
4034 #endif
4035                         if (erv <= -52)
4036                                 goto undfl;
4037                         LLval(&rv) = t00 >> (12 - erv);
4038                         Set_errno(ERANGE);
4039                         goto ret;
4040                         }
4041                 if (bexact) {
4042 #ifdef SET_INEXACT
4043                         if (!(t00 & 0x3ff) && !(tlo & 0xffffffffull)) {
4044                                 bc.inexact = 0;
4045                                 goto noround1;
4046                                 }
4047 #endif
4048 #ifdef Honor_FLT_ROUNDS
4049                         switch(bc.rounding) {
4050                           case 2:
4051                                 if (t00 & 0x3ff)
4052                                         goto roundup1;
4053                           case 0: goto noround1;
4054                           }
4055 #endif
4056                         if (t00 & 0x200 && (t00 & 0x5ff || tlo))
4057                                 goto roundup1;
4058                         goto noround1;
4059                         }
4060                 if ((tlo & 0xfffffff0) | (t00 & 0x1ff)
4061                  && (nd <= 19 ||  ((t00 + (1ull << i)) & 0x7ffffffffffffe00ull)
4062                                 == (t00 & 0x7ffffffffffffe00ull))) {
4063                         /* Unambiguous result. */
4064 #ifdef Honor_FLT_ROUNDS
4065                         switch(bc.rounding) {
4066                           case 0: goto noround1;
4067                           case 2: goto roundup1;
4068                           }
4069 #endif
4070                         if (t00 & 0x200) { /* round up */
4071  roundup1:
4072                                 t00 += 0x400;
4073                                 if (!(t00 & 0x4000000000000000ull)) {
4074                                         /* rounded up to a power of 2 */
4075                                         if (erv >= 0x7ff)
4076                                                 goto ovfl;
4077                                         terv = erv;
4078                                         LLval(&rv) = terv << 52;
4079                                         goto ret;
4080                                         }
4081                                 }
4082  noround1:
4083                         if (erv >= 0x800)
4084                                 goto ovfl;
4085                         terv = erv - 1;
4086                         LLval(&rv) = (terv << 52) | ((t00 & 0x3ffffffffffffc00ull) >> 10);
4087                         goto ret;
4088                         }
4089                 }
4090  many_digits:
4091         Debug(++dtoa_stats[2]);
4092         if (nd > 17) {
4093                 if (nd > 18) {
4094                         yz /= 100;
4095                         e1 += 2;
4096                         }
4097                 else {
4098                         yz /= 10;
4099                         e1 += 1;
4100                         }
4101                 y = yz / 100000000;
4102                 }
4103         else if (nd > 9) {
4104                 i = nd - 9;
4105                 y = (yz >> i) / pfive[i-1];
4106                 }
4107         else
4108                 y = yz;
4109         dval(&rv) = yz;
4110 #endif /*}*/
4111
4112 #ifdef IEEE_Arith
4113 #ifdef Avoid_Underflow
4114         bc.scale = 0;
4115 #endif
4116 #endif /*IEEE_Arith*/
4117
4118         /* Get starting approximation = rv * 10**e1 */
4119
4120         if (e1 > 0) {
4121                 if ((i = e1 & 15))
4122                         dval(&rv) *= tens[i];
4123                 if (e1 &= ~15) {
4124                         if (e1 > DBL_MAX_10_EXP) {
4125  ovfl:
4126                                 /* Can't trust HUGE_VAL */
4127 #ifdef IEEE_Arith
4128 #ifdef Honor_FLT_ROUNDS
4129                                 switch(bc.rounding) {
4130                                   case 0: /* toward 0 */
4131                                   case 3: /* toward -infinity */
4132                                         word0(&rv) = Big0;
4133                                         word1(&rv) = Big1;
4134                                         break;
4135                                   default:
4136                                         word0(&rv) = Exp_mask;
4137                                         word1(&rv) = 0;
4138                                   }
4139 #else /*Honor_FLT_ROUNDS*/
4140                                 word0(&rv) = Exp_mask;
4141                                 word1(&rv) = 0;
4142 #endif /*Honor_FLT_ROUNDS*/
4143 #ifdef SET_INEXACT
4144                                 /* set overflow bit */
4145                                 dval(&rv0) = 1e300;
4146                                 dval(&rv0) *= dval(&rv0);
4147 #endif
4148 #else /*IEEE_Arith*/
4149                                 word0(&rv) = Big0;
4150                                 word1(&rv) = Big1;
4151 #endif /*IEEE_Arith*/
4152  range_err:
4153                                 if (bd0) {
4154                                         Bfree(bb MTb);
4155                                         Bfree(bd MTb);
4156                                         Bfree(bs MTb);
4157                                         Bfree(bd0 MTb);
4158                                         Bfree(delta MTb);
4159                                         }
4160                                 Set_errno(ERANGE);
4161                                 goto ret;
4162                                 }
4163                         e1 >>= 4;
4164                         for(j = 0; e1 > 1; j++, e1 >>= 1)
4165                                 if (e1 & 1)
4166                                         dval(&rv) *= bigtens[j];
4167                 /* The last multiplication could overflow. */
4168                         word0(&rv) -= P*Exp_msk1;
4169                         dval(&rv) *= bigtens[j];
4170                         if ((z = word0(&rv) & Exp_mask)
4171                          > Exp_msk1*(DBL_MAX_EXP+Bias-P))
4172                                 goto ovfl;
4173                         if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
4174                                 /* set to largest number */
4175                                 /* (Can't trust DBL_MAX) */
4176                                 word0(&rv) = Big0;
4177                                 word1(&rv) = Big1;
4178                                 }
4179                         else
4180                                 word0(&rv) += P*Exp_msk1;
4181                         }
4182                 }
4183         else if (e1 < 0) {
4184                 e1 = -e1;
4185                 if ((i = e1 & 15))
4186                         dval(&rv) /= tens[i];
4187                 if (e1 >>= 4) {
4188                         if (e1 >= 1 << n_bigtens)
4189                                 goto undfl;
4190 #ifdef Avoid_Underflow
4191                         if (e1 & Scale_Bit)
4192                                 bc.scale = 2*P;
4193                         for(j = 0; e1 > 0; j++, e1 >>= 1)
4194                                 if (e1 & 1)
4195                                         dval(&rv) *= tinytens[j];
4196                         if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
4197                                                 >> Exp_shift)) > 0) {
4198                                 /* scaled rv is denormal; clear j low bits */
4199                                 if (j >= 32) {
4200                                         if (j > 54)
4201                                                 goto undfl;
4202                                         word1(&rv) = 0;
4203                                         if (j >= 53)
4204                                          word0(&rv) = (P+2)*Exp_msk1;
4205                                         else
4206                                          word0(&rv) &= 0xffffffff << (j-32);
4207                                         }
4208                                 else
4209                                         word1(&rv) &= 0xffffffff << j;
4210                                 }
4211 #else
4212                         for(j = 0; e1 > 1; j++, e1 >>= 1)
4213                                 if (e1 & 1)
4214                                         dval(&rv) *= tinytens[j];
4215                         /* The last multiplication could underflow. */
4216                         dval(&rv0) = dval(&rv);
4217                         dval(&rv) *= tinytens[j];
4218                         if (!dval(&rv)) {
4219                                 dval(&rv) = 2.*dval(&rv0);
4220                                 dval(&rv) *= tinytens[j];
4221 #endif
4222                                 if (!dval(&rv)) {
4223  undfl:
4224                                         dval(&rv) = 0.;
4225 #ifdef Honor_FLT_ROUNDS
4226                                         if (bc.rounding == 2)
4227                                                 word1(&rv) = 1;
4228 #endif
4229                                         goto range_err;
4230                                         }
4231 #ifndef Avoid_Underflow
4232                                 word0(&rv) = Tiny0;
4233                                 word1(&rv) = Tiny1;
4234                                 /* The refinement below will clean
4235                                  * this approximation up.
4236                                  */
4237                                 }
4238 #endif
4239                         }
4240                 }
4241
4242         /* Now the hard part -- adjusting rv to the correct value.*/
4243
4244         /* Put digits into bd: true value = bd * 10^e */
4245
4246         bc.nd = nd - nz1;
4247 #ifndef NO_STRTOD_BIGCOMP
4248         bc.nd0 = nd0;   /* Only needed if nd > strtod_diglim, but done here */
4249                         /* to silence an erroneous warning about bc.nd0 */
4250                         /* possibly not being initialized. */
4251         if (nd > strtod_diglim) {
4252                 /* ASSERT(strtod_diglim >= 18); 18 == one more than the */
4253                 /* minimum number of decimal digits to distinguish double values */
4254                 /* in IEEE arithmetic. */
4255                 i = j = 18;
4256                 if (i > nd0)
4257                         j += bc.dplen;
4258                 for(;;) {
4259                         if (--j < bc.dp1 && j >= bc.dp0)
4260                                 j = bc.dp0 - 1;
4261                         if (s0[j] != '0')
4262                                 break;
4263                         --i;
4264                         }
4265                 e += nd - i;
4266                 nd = i;
4267                 if (nd0 > nd)
4268                         nd0 = nd;
4269                 if (nd < 9) { /* must recompute y */
4270                         y = 0;
4271                         for(i = 0; i < nd0; ++i)
4272                                 y = 10*y + s0[i] - '0';
4273                         for(j = bc.dp1; i < nd; ++i)
4274                                 y = 10*y + s0[j++] - '0';
4275                         }
4276                 }
4277 #endif
4278         bd0 = s2b(s0, nd0, nd, y, bc.dplen MTb);
4279
4280         for(;;) {
4281                 bd = Balloc(bd0->k MTb);
4282                 Bcopy(bd, bd0);
4283                 bb = d2b(&rv, &bbe, &bbbits MTb);       /* rv = bb * 2^bbe */
4284                 bs = i2b(1 MTb);
4285
4286                 if (e >= 0) {
4287                         bb2 = bb5 = 0;
4288                         bd2 = bd5 = e;
4289                         }
4290                 else {
4291                         bb2 = bb5 = -e;
4292                         bd2 = bd5 = 0;
4293                         }
4294                 if (bbe >= 0)
4295                         bb2 += bbe;
4296                 else
4297                         bd2 -= bbe;
4298                 bs2 = bb2;
4299 #ifdef Honor_FLT_ROUNDS
4300                 if (bc.rounding != 1)
4301                         bs2++;
4302 #endif
4303 #ifdef Avoid_Underflow
4304                 Lsb = LSB;
4305                 Lsb1 = 0;
4306                 j = bbe - bc.scale;
4307                 i = j + bbbits - 1;     /* logb(rv) */
4308                 j = P + 1 - bbbits;
4309                 if (i < Emin) { /* denormal */
4310                         i = Emin - i;
4311                         j -= i;
4312                         if (i < 32)
4313                                 Lsb <<= i;
4314                         else if (i < 52)
4315                                 Lsb1 = Lsb << (i-32);
4316                         else
4317                                 Lsb1 = Exp_mask;
4318                         }
4319 #else /*Avoid_Underflow*/
4320 #ifdef Sudden_Underflow
4321 #ifdef IBM
4322                 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
4323 #else
4324                 j = P + 1 - bbbits;
4325 #endif
4326 #else /*Sudden_Underflow*/
4327                 j = bbe;
4328                 i = j + bbbits - 1;     /* logb(rv) */
4329                 if (i < Emin)   /* denormal */
4330                         j += P - Emin;
4331                 else
4332                         j = P + 1 - bbbits;
4333 #endif /*Sudden_Underflow*/
4334 #endif /*Avoid_Underflow*/
4335                 bb2 += j;
4336                 bd2 += j;
4337 #ifdef Avoid_Underflow
4338                 bd2 += bc.scale;
4339 #endif
4340                 i = bb2 < bd2 ? bb2 : bd2;
4341                 if (i > bs2)
4342                         i = bs2;
4343                 if (i > 0) {
4344                         bb2 -= i;
4345                         bd2 -= i;
4346                         bs2 -= i;
4347                         }
4348                 if (bb5 > 0) {
4349                         bs = pow5mult(bs, bb5 MTb);
4350                         bb1 = mult(bs, bb MTb);
4351                         Bfree(bb MTb);
4352                         bb = bb1;
4353                         }
4354                 if (bb2 > 0)
4355                         bb = lshift(bb, bb2 MTb);
4356                 if (bd5 > 0)
4357                         bd = pow5mult(bd, bd5 MTb);
4358                 if (bd2 > 0)
4359                         bd = lshift(bd, bd2 MTb);
4360                 if (bs2 > 0)
4361                         bs = lshift(bs, bs2 MTb);
4362                 delta = diff(bb, bd MTb);
4363                 bc.dsign = delta->sign;
4364                 delta->sign = 0;
4365                 i = cmp(delta, bs);
4366 #ifndef NO_STRTOD_BIGCOMP /*{*/
4367                 if (bc.nd > nd && i <= 0) {
4368                         if (bc.dsign) {
4369                                 /* Must use bigcomp(). */
4370                                 req_bigcomp = 1;
4371                                 break;
4372                                 }
4373 #ifdef Honor_FLT_ROUNDS
4374                         if (bc.rounding != 1) {
4375                                 if (i < 0) {
4376                                         req_bigcomp = 1;
4377                                         break;
4378                                         }
4379                                 }
4380                         else
4381 #endif
4382                                 i = -1; /* Discarded digits make delta smaller. */
4383                         }
4384 #endif /*}*/
4385 #ifdef Honor_FLT_ROUNDS /*{*/
4386                 if (bc.rounding != 1) {
4387                         if (i < 0) {
4388                                 /* Error is less than an ulp */
4389                                 if (!delta->x[0] && delta->wds <= 1) {
4390                                         /* exact */
4391 #ifdef SET_INEXACT
4392                                         bc.inexact = 0;
4393 #endif
4394                                         break;
4395                                         }
4396                                 if (bc.rounding) {
4397                                         if (bc.dsign) {
4398                                                 adj.d = 1.;
4399                                                 goto apply_adj;
4400                                                 }
4401                                         }
4402                                 else if (!bc.dsign) {
4403                                         adj.d = -1.;
4404                                         if (!word1(&rv)
4405                                          && !(word0(&rv) & Frac_mask)) {
4406                                                 y = word0(&rv) & Exp_mask;
4407 #ifdef Avoid_Underflow
4408                                                 if (!bc.scale || y > 2*P*Exp_msk1)
4409 #else
4410                                                 if (y)
4411 #endif
4412                                                   {
4413                                                   delta = lshift(delta,Log2P MTb);
4414                                                   if (cmp(delta, bs) <= 0)
4415                                                         adj.d = -0.5;
4416                                                   }
4417                                                 }
4418  apply_adj:
4419 #ifdef Avoid_Underflow /*{*/
4420                                         if (bc.scale && (y = word0(&rv) & Exp_mask)
4421                                                 <= 2*P*Exp_msk1)
4422                                           word0(&adj) += (2*P+1)*Exp_msk1 - y;
4423 #else
4424 #ifdef Sudden_Underflow
4425                                         if ((word0(&rv) & Exp_mask) <=
4426                                                         P*Exp_msk1) {
4427                                                 word0(&rv) += P*Exp_msk1;
4428                                                 dval(&rv) += adj.d*ulp(dval(&rv));
4429                                                 word0(&rv) -= P*Exp_msk1;
4430                                                 }
4431                                         else
4432 #endif /*Sudden_Underflow*/
4433 #endif /*Avoid_Underflow}*/
4434                                         dval(&rv) += adj.d*ulp(&rv);
4435                                         }
4436                                 break;
4437                                 }
4438                         adj.d = ratio(delta, bs);
4439                         if (adj.d < 1.)
4440                                 adj.d = 1.;
4441                         if (adj.d <= 0x7ffffffe) {
4442                                 /* adj = rounding ? ceil(adj) : floor(adj); */
4443                                 y = adj.d;
4444                                 if (y != adj.d) {
4445                                         if (!((bc.rounding>>1) ^ bc.dsign))
4446                                                 y++;
4447                                         adj.d = y;
4448                                         }
4449                                 }
4450 #ifdef Avoid_Underflow /*{*/
4451                         if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
4452                                 word0(&adj) += (2*P+1)*Exp_msk1 - y;
4453 #else
4454 #ifdef Sudden_Underflow
4455                         if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
4456                                 word0(&rv) += P*Exp_msk1;
4457                                 adj.d *= ulp(dval(&rv));
4458                                 if (bc.dsign)
4459                                         dval(&rv) += adj.d;
4460                                 else
4461                                         dval(&rv) -= adj.d;
4462                                 word0(&rv) -= P*Exp_msk1;
4463                                 goto cont;
4464                                 }
4465 #endif /*Sudden_Underflow*/
4466 #endif /*Avoid_Underflow}*/
4467                         adj.d *= ulp(&rv);
4468                         if (bc.dsign) {
4469                                 if (word0(&rv) == Big0 && word1(&rv) == Big1)
4470                                         goto ovfl;
4471                                 dval(&rv) += adj.d;
4472                                 }
4473                         else
4474                                 dval(&rv) -= adj.d;
4475                         goto cont;
4476                         }
4477 #endif /*}Honor_FLT_ROUNDS*/
4478
4479                 if (i < 0) {
4480                         /* Error is less than half an ulp -- check for
4481                          * special case of mantissa a power of two.
4482                          */
4483                         if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask
4484 #ifdef IEEE_Arith /*{*/
4485 #ifdef Avoid_Underflow
4486                          || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
4487 #else
4488                          || (word0(&rv) & Exp_mask) <= Exp_msk1
4489 #endif
4490 #endif /*}*/
4491                                 ) {
4492 #ifdef SET_INEXACT
4493                                 if (!delta->x[0] && delta->wds <= 1)
4494                                         bc.inexact = 0;
4495 #endif
4496                                 break;
4497                                 }
4498                         if (!delta->x[0] && delta->wds <= 1) {
4499                                 /* exact result */
4500 #ifdef SET_INEXACT
4501                                 bc.inexact = 0;
4502 #endif
4503                                 break;
4504                                 }
4505                         delta = lshift(delta,Log2P MTb);
4506                         if (cmp(delta, bs) > 0)
4507                                 goto drop_down;
4508                         break;
4509                         }
4510                 if (i == 0) {
4511                         /* exactly half-way between */
4512                         if (bc.dsign) {
4513                                 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
4514                                  &&  word1(&rv) == (
4515 #ifdef Avoid_Underflow
4516                         (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
4517                 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
4518 #endif
4519                                                    0xffffffff)) {
4520                                         /*boundary case -- increment exponent*/
4521                                         if (word0(&rv) == Big0 && word1(&rv) == Big1)
4522                                                 goto ovfl;
4523                                         word0(&rv) = (word0(&rv) & Exp_mask)
4524                                                 + Exp_msk1
4525 #ifdef IBM
4526                                                 | Exp_msk1 >> 4
4527 #endif
4528                                                 ;
4529                                         word1(&rv) = 0;
4530 #ifdef Avoid_Underflow
4531                                         bc.dsign = 0;
4532 #endif
4533                                         break;
4534                                         }
4535                                 }
4536                         else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
4537  drop_down:
4538                                 /* boundary case -- decrement exponent */
4539 #ifdef Sudden_Underflow /*{{*/
4540                                 L = word0(&rv) & Exp_mask;
4541 #ifdef IBM
4542                                 if (L <  Exp_msk1)
4543 #else
4544 #ifdef Avoid_Underflow
4545                                 if (L <= (bc.scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
4546 #else
4547                                 if (L <= Exp_msk1)
4548 #endif /*Avoid_Underflow*/
4549 #endif /*IBM*/
4550                                         {
4551                                         if (bc.nd >nd) {
4552                                                 bc.uflchk = 1;
4553                                                 break;
4554                                                 }
4555                                         goto undfl;
4556                                         }
4557                                 L -= Exp_msk1;
4558 #else /*Sudden_Underflow}{*/
4559 #ifdef Avoid_Underflow
4560                                 if (bc.scale) {
4561                                         L = word0(&rv) & Exp_mask;
4562                                         if (L <= (2*P+1)*Exp_msk1) {
4563                                                 if (L > (P+2)*Exp_msk1)
4564                                                         /* round even ==> */
4565                                                         /* accept rv */
4566                                                         break;
4567                                                 /* rv = smallest denormal */
4568                                                 if (bc.nd >nd) {
4569                                                         bc.uflchk = 1;
4570                                                         break;
4571                                                         }
4572                                                 goto undfl;
4573                                                 }
4574                                         }
4575 #endif /*Avoid_Underflow*/
4576                                 L = (word0(&rv) & Exp_mask) - Exp_msk1;
4577 #endif /*Sudden_Underflow}}*/
4578                                 word0(&rv) = L | Bndry_mask1;
4579                                 word1(&rv) = 0xffffffff;
4580 #ifdef IBM
4581                                 goto cont;
4582 #else
4583 #ifndef NO_STRTOD_BIGCOMP
4584                                 if (bc.nd > nd)
4585                                         goto cont;
4586 #endif
4587                                 break;
4588 #endif
4589                                 }
4590 #ifndef ROUND_BIASED
4591 #ifdef Avoid_Underflow
4592                         if (Lsb1) {
4593                                 if (!(word0(&rv) & Lsb1))
4594                                         break;
4595                                 }
4596                         else if (!(word1(&rv) & Lsb))
4597                                 break;
4598 #else
4599                         if (!(word1(&rv) & LSB))
4600                                 break;
4601 #endif
4602 #endif
4603                         if (bc.dsign)
4604 #ifdef Avoid_Underflow
4605                                 dval(&rv) += sulp(&rv, &bc);
4606 #else
4607                                 dval(&rv) += ulp(&rv);
4608 #endif
4609 #ifndef ROUND_BIASED
4610                         else {
4611 #ifdef Avoid_Underflow
4612                                 dval(&rv) -= sulp(&rv, &bc);
4613 #else
4614                                 dval(&rv) -= ulp(&rv);
4615 #endif
4616 #ifndef Sudden_Underflow
4617                                 if (!dval(&rv)) {
4618                                         if (bc.nd >nd) {
4619                                                 bc.uflchk = 1;
4620                                                 break;
4621                                                 }
4622                                         goto undfl;
4623                                         }
4624 #endif
4625                                 }
4626 #ifdef Avoid_Underflow
4627                         bc.dsign = 1 - bc.dsign;
4628 #endif
4629 #endif
4630                         break;
4631                         }
4632                 if ((aadj = ratio(delta, bs)) <= 2.) {
4633                         if (bc.dsign)
4634                                 aadj = aadj1 = 1.;
4635                         else if (word1(&rv) || word0(&rv) & Bndry_mask) {
4636 #ifndef Sudden_Underflow
4637                                 if (word1(&rv) == Tiny1 && !word0(&rv)) {
4638                                         if (bc.nd >nd) {
4639                                                 bc.uflchk = 1;
4640                                                 break;
4641                                                 }
4642                                         goto undfl;
4643                                         }
4644 #endif
4645                                 aadj = 1.;
4646                                 aadj1 = -1.;
4647                                 }
4648                         else {
4649                                 /* special case -- power of FLT_RADIX to be */
4650                                 /* rounded down... */
4651
4652                                 if (aadj < 2./FLT_RADIX)
4653                                         aadj = 1./FLT_RADIX;
4654                                 else
4655                                         aadj *= 0.5;
4656                                 aadj1 = -aadj;
4657                                 }
4658                         }
4659                 else {
4660                         aadj *= 0.5;
4661                         aadj1 = bc.dsign ? aadj : -aadj;
4662 #ifdef Check_FLT_ROUNDS
4663                         switch(bc.rounding) {
4664                                 case 2: /* towards +infinity */
4665                                         aadj1 -= 0.5;
4666                                         break;
4667                                 case 0: /* towards 0 */
4668                                 case 3: /* towards -infinity */
4669                                         aadj1 += 0.5;
4670                                 }
4671 #else
4672                         if (Flt_Rounds == 0)
4673                                 aadj1 += 0.5;
4674 #endif /*Check_FLT_ROUNDS*/
4675                         }
4676                 y = word0(&rv) & Exp_mask;
4677
4678                 /* Check for overflow */
4679
4680                 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
4681                         dval(&rv0) = dval(&rv);
4682                         word0(&rv) -= P*Exp_msk1;
4683                         adj.d = aadj1 * ulp(&rv);
4684                         dval(&rv) += adj.d;
4685                         if ((word0(&rv) & Exp_mask) >=
4686                                         Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
4687                                 if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
4688                                         goto ovfl;
4689                                 word0(&rv) = Big0;
4690                                 word1(&rv) = Big1;
4691                                 goto cont;
4692                                 }
4693                         else
4694                                 word0(&rv) += P*Exp_msk1;
4695                         }
4696                 else {
4697 #ifdef Avoid_Underflow
4698                         if (bc.scale && y <= 2*P*Exp_msk1) {
4699                                 if (aadj <= 0x7fffffff) {
4700                                         if ((z = aadj) <= 0)
4701                                                 z = 1;
4702                                         aadj = z;
4703                                         aadj1 = bc.dsign ? aadj : -aadj;
4704                                         }
4705                                 dval(&aadj2) = aadj1;
4706                                 word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
4707                                 aadj1 = dval(&aadj2);
4708                                 adj.d = aadj1 * ulp(&rv);
4709                                 dval(&rv) += adj.d;
4710                                 if (rv.d == 0.)
4711 #ifdef NO_STRTOD_BIGCOMP
4712                                         goto undfl;
4713 #else
4714                                         {
4715                                         req_bigcomp = 1;
4716                                         break;
4717                                         }
4718 #endif
4719                                 }
4720                         else {
4721                                 adj.d = aadj1 * ulp(&rv);
4722                                 dval(&rv) += adj.d;
4723                                 }
4724 #else
4725 #ifdef Sudden_Underflow
4726                         if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
4727                                 dval(&rv0) = dval(&rv);
4728                                 word0(&rv) += P*Exp_msk1;
4729                                 adj.d = aadj1 * ulp(&rv);
4730                                 dval(&rv) += adj.d;
4731 #ifdef IBM
4732                                 if ((word0(&rv) & Exp_mask) <  P*Exp_msk1)
4733 #else
4734                                 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
4735 #endif
4736                                         {
4737                                         if (word0(&rv0) == Tiny0
4738                                          && word1(&rv0) == Tiny1) {
4739                                                 if (bc.nd >nd) {
4740                                                         bc.uflchk = 1;
4741                                                         break;
4742                                                         }
4743                                                 goto undfl;
4744                                                 }
4745                                         word0(&rv) = Tiny0;
4746                                         word1(&rv) = Tiny1;
4747                                         goto cont;
4748                                         }
4749                                 else
4750                                         word0(&rv) -= P*Exp_msk1;
4751                                 }
4752                         else {
4753                                 adj.d = aadj1 * ulp(&rv);
4754                                 dval(&rv) += adj.d;
4755                                 }
4756 #else /*Sudden_Underflow*/
4757                         /* Compute adj so that the IEEE rounding rules will
4758                          * correctly round rv + adj in some half-way cases.
4759                          * If rv * ulp(rv) is denormalized (i.e.,
4760                          * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
4761                          * trouble from bits lost to denormalization;
4762                          * example: 1.2e-307 .
4763                          */
4764                         if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
4765                                 aadj1 = (double)(int)(aadj + 0.5);
4766                                 if (!bc.dsign)
4767                                         aadj1 = -aadj1;
4768                                 }
4769                         adj.d = aadj1 * ulp(&rv);
4770                         dval(&rv) += adj.d;
4771 #endif /*Sudden_Underflow*/
4772 #endif /*Avoid_Underflow*/
4773                         }
4774                 z = word0(&rv) & Exp_mask;
4775 #ifndef SET_INEXACT
4776                 if (bc.nd == nd) {
4777 #ifdef Avoid_Underflow
4778                 if (!bc.scale)
4779 #endif
4780                 if (y == z) {
4781                         /* Can we stop now? */
4782                         L = (Long)aadj;
4783                         aadj -= L;
4784                         /* The tolerances below are conservative. */
4785                         if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
4786                                 if (aadj < .4999999 || aadj > .5000001)
4787                                         break;
4788                                 }
4789                         else if (aadj < .4999999/FLT_RADIX)
4790                                 break;
4791                         }
4792                 }
4793 #endif
4794  cont:
4795                 Bfree(bb MTb);
4796                 Bfree(bd MTb);
4797                 Bfree(bs MTb);
4798                 Bfree(delta MTb);
4799                 }
4800         Bfree(bb MTb);
4801         Bfree(bd MTb);
4802         Bfree(bs MTb);
4803         Bfree(bd0 MTb);
4804         Bfree(delta MTb);
4805 #ifndef NO_STRTOD_BIGCOMP
4806         if (req_bigcomp) {
4807                 bd0 = 0;
4808                 bc.e0 += nz1;
4809                 bigcomp(&rv, s0, &bc MTb);
4810                 y = word0(&rv) & Exp_mask;
4811                 if (y == Exp_mask)
4812                         goto ovfl;
4813                 if (y == 0 && rv.d == 0.)
4814                         goto undfl;
4815                 }
4816 #endif
4817 #ifdef Avoid_Underflow
4818         if (bc.scale) {
4819                 word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
4820                 word1(&rv0) = 0;
4821                 dval(&rv) *= dval(&rv0);
4822 #ifndef NO_ERRNO
4823                 /* try to avoid the bug of testing an 8087 register value */
4824 #ifdef IEEE_Arith
4825                 if (!(word0(&rv) & Exp_mask))
4826 #else
4827                 if (word0(&rv) == 0 && word1(&rv) == 0)
4828 #endif
4829                         Set_errno(ERANGE);
4830 #endif
4831                 }
4832 #endif /* Avoid_Underflow */
4833  ret:
4834 #ifdef SET_INEXACT
4835         if (bc.inexact) {
4836                 if (!(word0(&rv) & Exp_mask)) {
4837                         /* set underflow and inexact bits */
4838                         dval(&rv0) = 1e-300;
4839                         dval(&rv0) *= dval(&rv0);
4840                         }
4841                 else if (!oldinexact) {
4842                         word0(&rv0) = Exp_1 + (70 << Exp_shift);
4843                         word1(&rv0) = 0;
4844                         dval(&rv0) += 1.;
4845                         }
4846                 }
4847         else if (!oldinexact)
4848                 clear_inexact();
4849 #endif
4850         if (se)
4851                 *se = (char *)s;
4852         return sign ? -dval(&rv) : dval(&rv);
4853         }
4854
4855 #ifndef MULTIPLE_THREADS
4856  static char *dtoa_result;
4857 #endif
4858
4859  static char *
4860 rv_alloc(int i MTd)
4861 {
4862         int j, k, *r;
4863
4864         j = sizeof(ULong);
4865         for(k = 0;
4866                 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
4867                 j <<= 1)
4868                         k++;
4869         r = (int*)Balloc(k MTa);
4870         *r = k;
4871         return
4872 #ifndef MULTIPLE_THREADS
4873         dtoa_result =
4874 #endif
4875                 (char *)(r+1);
4876         }
4877
4878  static char *
4879 nrv_alloc(const char *s, char *s0, size_t s0len, char **rve, int n MTd)
4880 {
4881         char *rv, *t;
4882
4883         if (!s0)
4884                 s0 = rv_alloc(n MTa);
4885         else if (s0len <= n) {
4886                 rv = 0;
4887                 t = rv + n;
4888                 goto rve_chk;
4889                 }
4890         t = rv = s0;
4891         while((*t = *s++))
4892                 ++t;
4893  rve_chk:
4894         if (rve)
4895                 *rve = t;
4896         return rv;
4897         }
4898
4899 /* freedtoa(s) must be used to free values s returned by dtoa
4900  * when MULTIPLE_THREADS is #defined.  It should be used in all cases,
4901  * but for consistency with earlier versions of dtoa, it is optional
4902  * when MULTIPLE_THREADS is not defined.
4903  */
4904
4905  void
4906 freedtoa(char *s)
4907 {
4908 #ifdef MULTIPLE_THREADS
4909         ThInfo *TI = 0;
4910 #endif
4911         Bigint *b = (Bigint *)((int *)s - 1);
4912         b->maxwds = 1 << (b->k = *(int*)b);
4913         Bfree(b MTb);
4914 #ifndef MULTIPLE_THREADS
4915         if (s == dtoa_result)
4916                 dtoa_result = 0;
4917 #endif
4918         }
4919
4920 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
4921  *
4922  * Inspired by "How to Print Floating-Point Numbers Accurately" by
4923  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
4924  *
4925  * Modifications:
4926  *      1. Rather than iterating, we use a simple numeric overestimate
4927  *         to determine k = floor(log10(d)).  We scale relevant
4928  *         quantities using O(log2(k)) rather than O(k) multiplications.
4929  *      2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
4930  *         try to generate digits strictly left to right.  Instead, we
4931  *         compute with fewer bits and propagate the carry if necessary
4932  *         when rounding the final digit up.  This is often faster.
4933  *      3. Under the assumption that input will be rounded nearest,
4934  *         mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
4935  *         That is, we allow equality in stopping tests when the
4936  *         round-nearest rule will give the same floating-point value
4937  *         as would satisfaction of the stopping test with strict
4938  *         inequality.
4939  *      4. We remove common factors of powers of 2 from relevant
4940  *         quantities.
4941  *      5. When converting floating-point integers less than 1e16,
4942  *         we use floating-point arithmetic rather than resorting
4943  *         to multiple-precision integers.
4944  *      6. When asked to produce fewer than 15 digits, we first try
4945  *         to get by with floating-point arithmetic; we resort to
4946  *         multiple-precision integer arithmetic only if we cannot
4947  *         guarantee that the floating-point calculation has given
4948  *         the correctly rounded result.  For k requested digits and
4949  *         "uniformly" distributed input, the probability is
4950  *         something like 10^(k-15) that we must resort to the Long
4951  *         calculation.
4952  */
4953
4954  char *
4955 dtoa_r(double dd, int mode, int ndigits, int *decpt, int *sign, char **rve, char *buf, size_t blen)
4956 {
4957  /*     Arguments ndigits, decpt, sign are similar to those
4958         of ecvt and fcvt; trailing zeros are suppressed from
4959         the returned string.  If not null, *rve is set to point
4960         to the end of the return value.  If d is +-Infinity or NaN,
4961         then *decpt is set to 9999.
4962
4963         mode:
4964                 0 ==> shortest string that yields d when read in
4965                         and rounded to nearest.
4966                 1 ==> like 0, but with Steele & White stopping rule;
4967                         e.g. with IEEE P754 arithmetic , mode 0 gives
4968                         1e23 whereas mode 1 gives 9.999999999999999e22.
4969                 2 ==> max(1,ndigits) significant digits.  This gives a
4970                         return value similar to that of ecvt, except
4971                         that trailing zeros are suppressed.
4972                 3 ==> through ndigits past the decimal point.  This
4973                         gives a return value similar to that from fcvt,
4974                         except that trailing zeros are suppressed, and
4975                         ndigits can be negative.
4976                 4,5 ==> similar to 2 and 3, respectively, but (in
4977                         round-nearest mode) with the tests of mode 0 to
4978                         possibly return a shorter string that rounds to d.
4979                         With IEEE arithmetic and compilation with
4980                         -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
4981                         as modes 2 and 3 when FLT_ROUNDS != 1.
4982                 6-9 ==> Debugging modes similar to mode - 4:  don't try
4983                         fast floating-point estimate (if applicable).
4984
4985                 Values of mode other than 0-9 are treated as mode 0.
4986
4987         When not NULL, buf is an output buffer of length blen, which must
4988         be large enough to accommodate suppressed trailing zeros and a trailing
4989         null byte.  If blen is too small, rv = NULL is returned, in which case
4990         if rve is not NULL, a subsequent call with blen >= (*rve - rv) + 1
4991         should succeed in returning buf.
4992
4993         When buf is NULL, sufficient space is allocated for the return value,
4994         which, when done using, the caller should pass to freedtoa().
4995
4996         USE_BF is automatically defined when neither NO_LONG_LONG nor NO_BF96
4997         is defined.
4998         */
4999
5000 #ifdef MULTIPLE_THREADS
5001         ThInfo *TI = 0;
5002 #endif
5003         int bbits, b2, b5, be, dig, i, ilim, ilim1,
5004                 j, j1, k, leftright, m2, m5, s2, s5, spec_case;
5005 #if !defined(Sudden_Underflow) || defined(USE_BF96)
5006         int denorm;
5007 #endif
5008         Bigint *b, *b1, *delta, *mlo, *mhi, *S;
5009         U u;
5010         char *s;
5011 #ifdef SET_INEXACT
5012         int inexact, oldinexact;
5013 #endif
5014 #ifdef USE_BF96 /*{{*/
5015         BF96 *p10;
5016         ULLong dbhi, dbits, dblo, den, hb, rb, rblo, res, res0, res3, reslo, sres,
5017                 sulp, tv0, tv1, tv2, tv3, ulp, ulplo, ulpmask, ures, ureslo, zb;
5018         int eulp, k1, n2, ulpadj, ulpshift;
5019 #else /*}{*/
5020 #ifndef Sudden_Underflow
5021         ULong x;
5022 #endif
5023         Long L;
5024         U d2, eps;
5025         double ds;
5026         int ieps, ilim0, k0, k_check, try_quick;
5027 #ifndef No_leftright
5028 #ifdef IEEE_Arith
5029         U eps1;
5030 #endif
5031 #endif
5032 #endif /*}}*/
5033 #ifdef Honor_FLT_ROUNDS /*{*/
5034         int Rounding;
5035 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
5036         Rounding = Flt_Rounds;
5037 #else /*}{*/
5038         Rounding = 1;
5039         switch(fegetround()) {
5040           case FE_TOWARDZERO:   Rounding = 0; break;
5041           case FE_UPWARD:       Rounding = 2; break;
5042           case FE_DOWNWARD:     Rounding = 3;
5043           }
5044 #endif /*}}*/
5045 #endif /*}*/
5046
5047         u.d = dd;
5048         if (word0(&u) & Sign_bit) {
5049                 /* set sign for everything, including 0's and NaNs */
5050                 *sign = 1;
5051                 word0(&u) &= ~Sign_bit; /* clear sign bit */
5052                 }
5053         else
5054                 *sign = 0;
5055
5056 #if defined(IEEE_Arith) + defined(VAX)
5057 #ifdef IEEE_Arith
5058         if ((word0(&u) & Exp_mask) == Exp_mask)
5059 #else
5060         if (word0(&u)  == 0x8000)
5061 #endif
5062                 {
5063                 /* Infinity or NaN */
5064                 *decpt = 9999;
5065 #ifdef IEEE_Arith
5066                 if (!word1(&u) && !(word0(&u) & 0xfffff))
5067                         return nrv_alloc("Infinity", buf, blen, rve, 8 MTb);
5068 #endif
5069                 return nrv_alloc("NaN", buf, blen, rve, 3 MTb);
5070                 }
5071 #endif
5072 #ifdef IBM
5073         dval(&u) += 0; /* normalize */
5074 #endif
5075         if (!dval(&u)) {
5076                 *decpt = 1;
5077                 return nrv_alloc("0", buf, blen, rve, 1 MTb);
5078                 }
5079
5080 #ifdef SET_INEXACT
5081 #ifndef USE_BF96
5082         try_quick =
5083 #endif
5084         oldinexact = get_inexact();
5085         inexact = 1;
5086 #endif
5087 #ifdef Honor_FLT_ROUNDS
5088         if (Rounding >= 2) {
5089                 if (*sign)
5090                         Rounding = Rounding == 2 ? 0 : 2;
5091                 else
5092                         if (Rounding != 2)
5093                                 Rounding = 0;
5094                 }
5095 #endif
5096 #ifdef USE_BF96 /*{{*/
5097         dbits = (u.LL & 0xfffffffffffffull) << 11;      /* fraction bits */
5098         if ((be = u.LL >> 52)) /* biased exponent; nonzero ==> normal */ {
5099                 dbits |= 0x8000000000000000ull;
5100                 denorm = ulpadj = 0;
5101                 }
5102         else {
5103                 denorm = 1;
5104                 ulpadj = be + 1;
5105                 dbits <<= 1;
5106                 if (!(dbits & 0xffffffff00000000ull)) {
5107                         dbits <<= 32;
5108                         be -= 32;
5109                         }
5110                 if (!(dbits & 0xffff000000000000ull)) {
5111                         dbits <<= 16;
5112                         be -= 16;
5113                         }
5114                 if (!(dbits & 0xff00000000000000ull)) {
5115                         dbits <<= 8;
5116                         be -= 8;
5117                         }
5118                 if (!(dbits & 0xf000000000000000ull)) {
5119                         dbits <<= 4;
5120                         be -= 4;
5121                         }
5122                 if (!(dbits & 0xc000000000000000ull)) {
5123                         dbits <<= 2;
5124                         be -= 2;
5125                         }
5126                 if (!(dbits & 0x8000000000000000ull)) {
5127                         dbits <<= 1;
5128                         be -= 1;
5129                         }
5130                 assert(be >= -51);
5131                 ulpadj -= be;
5132                 }
5133         j = Lhint[be + 51];
5134         p10 = &pten[j];
5135         dbhi = dbits >> 32;
5136         dblo = dbits & 0xffffffffull;
5137         i = be - 0x3fe;
5138         if (i < p10->e
5139         || (i == p10->e && (dbhi < p10->b0 || (dbhi == p10->b0 && dblo < p10->b1))))
5140                 --j;
5141         k = j - 342;
5142
5143         /* now 10^k <= dd < 10^(k+1) */
5144
5145 #else /*}{*/
5146
5147         b = d2b(&u, &be, &bbits MTb);
5148 #ifdef Sudden_Underflow
5149         i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
5150 #else
5151         if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
5152 #endif
5153                 dval(&d2) = dval(&u);
5154                 word0(&d2) &= Frac_mask1;
5155                 word0(&d2) |= Exp_11;
5156 #ifdef IBM
5157                 if (j = 11 - hi0bits(word0(&d2) & Frac_mask))
5158                         dval(&d2) /= 1 << j;
5159 #endif
5160
5161                 /* log(x)       ~=~ log(1.5) + (x-1.5)/1.5
5162                  * log10(x)      =  log(x) / log(10)
5163                  *              ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
5164                  * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
5165                  *
5166                  * This suggests computing an approximation k to log10(d) by
5167                  *
5168                  * k = (i - Bias)*0.301029995663981
5169                  *      + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
5170                  *
5171                  * We want k to be too large rather than too small.
5172                  * The error in the first-order Taylor series approximation
5173                  * is in our favor, so we just round up the constant enough
5174                  * to compensate for any error in the multiplication of
5175                  * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
5176                  * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
5177                  * adding 1e-13 to the constant term more than suffices.
5178                  * Hence we adjust the constant term to 0.1760912590558.
5179                  * (We could get a more accurate k by invoking log10,
5180                  *  but this is probably not worthwhile.)
5181                  */
5182
5183                 i -= Bias;
5184 #ifdef IBM
5185                 i <<= 2;
5186                 i += j;
5187 #endif
5188 #ifndef Sudden_Underflow
5189                 denorm = 0;
5190                 }
5191         else {
5192                 /* d is denormalized */
5193
5194                 i = bbits + be + (Bias + (P-1) - 1);
5195                 x = i > 32  ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
5196                             : word1(&u) << (32 - i);
5197                 dval(&d2) = x;
5198                 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
5199                 i -= (Bias + (P-1) - 1) + 1;
5200                 denorm = 1;
5201                 }
5202 #endif
5203         ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
5204         k = (int)ds;
5205         if (ds < 0. && ds != k)
5206                 k--;    /* want k = floor(ds) */
5207         k_check = 1;
5208         if (k >= 0 && k <= Ten_pmax) {
5209                 if (dval(&u) < tens[k])
5210                         k--;
5211                 k_check = 0;
5212                 }
5213         j = bbits - i - 1;
5214         if (j >= 0) {
5215                 b2 = 0;
5216                 s2 = j;
5217                 }
5218         else {
5219                 b2 = -j;
5220                 s2 = 0;
5221                 }
5222         if (k >= 0) {
5223                 b5 = 0;
5224                 s5 = k;
5225                 s2 += k;
5226                 }
5227         else {
5228                 b2 -= k;
5229                 b5 = -k;
5230                 s5 = 0;
5231                 }
5232 #endif /*}}*/
5233         if (mode < 0 || mode > 9)
5234                 mode = 0;
5235
5236 #ifndef USE_BF96
5237 #ifndef SET_INEXACT
5238 #ifdef Check_FLT_ROUNDS
5239         try_quick = Rounding == 1;
5240 #endif
5241 #endif /*SET_INEXACT*/
5242 #endif
5243
5244         if (mode > 5) {
5245                 mode -= 4;
5246 #ifndef USE_BF96
5247                 try_quick = 0;
5248 #endif
5249                 }
5250         leftright = 1;
5251         ilim = ilim1 = -1;      /* Values for cases 0 and 1; done here to */
5252                                 /* silence erroneous "gcc -Wall" warning. */
5253         switch(mode) {
5254                 case 0:
5255                 case 1:
5256                         i = 18;
5257                         ndigits = 0;
5258                         break;
5259                 case 2:
5260                         leftright = 0;
5261                         /* no break */
5262                 case 4:
5263                         if (ndigits <= 0)
5264                                 ndigits = 1;
5265                         ilim = ilim1 = i = ndigits;
5266                         break;
5267                 case 3:
5268                         leftright = 0;
5269                         /* no break */
5270                 case 5:
5271                         i = ndigits + k + 1;
5272                         ilim = i;
5273                         ilim1 = i - 1;
5274                         if (i <= 0)
5275                                 i = 1;
5276                 }
5277         if (!buf) {
5278        &n