0029399: Optimize reading of floating point values from text strings -- base dtoa.c
[occt.git] / src / Standard / Standard_Strtod.cxx
1 /****************************************************************
2  *
3  * The author of this software is David M. Gay.
4  *
5  * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
6  *
7  * Permission to use, copy, modify, and distribute this software for any
8  * purpose without fee is hereby granted, provided that this entire notice
9  * is included in all copies of any software which is or includes a copy
10  * or modification of this software and in all copies of the supporting
11  * documentation for such software.
12  *
13  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
14  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
15  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
16  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
17  *
18  ***************************************************************/
19
20 /* Please send bug reports to David M. Gay (dmg at acm dot org,
21  * with " at " changed at "@" and " dot " changed to ".").      */
22
23 /* On a machine with IEEE extended-precision registers, it is
24  * necessary to specify double-precision (53-bit) rounding precision
25  * before invoking strtod or dtoa.  If the machine uses (the equivalent
26  * of) Intel 80x87 arithmetic, the call
27  *      _control87(PC_53, MCW_PC);
28  * does this with many compilers.  Whether this or another call is
29  * appropriate depends on the compiler; for this to work, it may be
30  * necessary to #include "float.h" or another system-dependent header
31  * file.
32  */
33
34 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
35  * (Note that IEEE arithmetic is disabled by gcc's -ffast-math flag.)
36  *
37  * This strtod returns a nearest machine number to the input decimal
38  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
39  * broken by the IEEE round-even rule.  Otherwise ties are broken by
40  * biased rounding (add half and chop).
41  *
42  * Inspired loosely by William D. Clinger's paper "How to Read Floating
43  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
44  *
45  * Modifications:
46  *
47  *      1. We only require IEEE, IBM, or VAX double-precision
48  *              arithmetic (not IEEE double-extended).
49  *      2. We get by with floating-point arithmetic in a case that
50  *              Clinger missed -- when we're computing d * 10^n
51  *              for a small integer d and the integer n is not too
52  *              much larger than 22 (the maximum integer k for which
53  *              we can represent 10^k exactly), we may be able to
54  *              compute (d*10^k) * 10^(e-k) with just one roundoff.
55  *      3. Rather than a bit-at-a-time adjustment of the binary
56  *              result in the hard case, we use floating-point
57  *              arithmetic to determine the adjustment to within
58  *              one bit; only in really hard cases do we need to
59  *              compute a second residual.
60  *      4. Because of 3., we don't need a large table of powers of 10
61  *              for ten-to-e (just some small tables, e.g. of 10^k
62  *              for 0 <= k <= 22).
63  */
64
65 /*
66  * #define IEEE_8087 for IEEE-arithmetic machines where the least
67  *      significant byte has the lowest address.
68  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
69  *      significant byte has the lowest address.
70  * #define Long int on machines with 32-bit ints and 64-bit longs.
71  * #define IBM for IBM mainframe-style floating-point arithmetic.
72  * #define VAX for VAX-style floating-point arithmetic (D_floating).
73  * #define No_leftright to omit left-right logic in fast floating-point
74  *      computation of dtoa.  This will cause dtoa modes 4 and 5 to be
75  *      treated the same as modes 2 and 3 for some inputs.
76  * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
77  *      and strtod and dtoa should round accordingly.  Unless Trust_FLT_ROUNDS
78  *      is also #defined, fegetround() will be queried for the rounding mode.
79  *      Note that both FLT_ROUNDS and fegetround() are specified by the C99
80  *      standard (and are specified to be consistent, with fesetround()
81  *      affecting the value of FLT_ROUNDS), but that some (Linux) systems
82  *      do not work correctly in this regard, so using fegetround() is more
83  *      portable than using FLT_ROUNDS directly.
84  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
85  *      and Honor_FLT_ROUNDS is not #defined.
86  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
87  *      that use extended-precision instructions to compute rounded
88  *      products and quotients) with IBM.
89  * #define ROUND_BIASED for IEEE-format with biased rounding and arithmetic
90  *      that rounds toward +Infinity.
91  * #define ROUND_BIASED_without_Round_Up for IEEE-format with biased
92  *      rounding when the underlying floating-point arithmetic uses
93  *      unbiased rounding.  This prevent using ordinary floating-point
94  *      arithmetic when the result could be computed with one rounding error.
95  * #define Inaccurate_Divide for IEEE-format with correctly rounded
96  *      products but inaccurate quotients, e.g., for Intel i860.
97  * #define NO_LONG_LONG on machines that do not have a "long long"
98  *      integer type (of >= 64 bits).  On such machines, you can
99  *      #define Just_16 to store 16 bits per 32-bit Long when doing
100  *      high-precision integer arithmetic.  Whether this speeds things
101  *      up or slows things down depends on the machine and the number
102  *      being converted.  If long long is available and the name is
103  *      something other than "long long", #define Llong to be the name,
104  *      and if "unsigned Llong" does not work as an unsigned version of
105  *      Llong, #define #ULLong to be the corresponding unsigned type.
106  * #define Bad_float_h if your system lacks a float.h or if it does not
107  *      define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
108  *      FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
109  * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
110  *      if memory is available and otherwise does something you deem
111  *      appropriate.  If MALLOC is undefined, malloc will be invoked
112  *      directly -- and assumed always to succeed.  Similarly, if you
113  *      want something other than the system's free() to be called to
114  *      recycle memory acquired from MALLOC, #define FREE to be the
115  *      name of the alternate routine.  (FREE or free is only called in
116  *      pathological cases, e.g., in a dtoa call after a dtoa return in
117  *      mode 3 with thousands of digits requested.)
118  * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
119  *      memory allocations from a private pool of memory when possible.
120  *      When used, the private pool is PRIVATE_MEM bytes long:  2304 bytes,
121  *      unless #defined to be a different length.  This default length
122  *      suffices to get rid of MALLOC calls except for unusual cases,
123  *      such as decimal-to-binary conversion of a very long string of
124  *      digits.  The longest string dtoa can return is about 751 bytes
125  *      long.  For conversions by strtod of strings of 800 digits and
126  *      all dtoa conversions in single-threaded executions with 8-byte
127  *      pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
128  *      pointers, PRIVATE_MEM >= 7112 appears adequate.
129  * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
130  *      #defined automatically on IEEE systems.  On such systems,
131  *      when INFNAN_CHECK is #defined, strtod checks
132  *      for Infinity and NaN (case insensitively).  On some systems
133  *      (e.g., some HP systems), it may be necessary to #define NAN_WORD0
134  *      appropriately -- to the most significant word of a quiet NaN.
135  *      (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
136  *      When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
137  *      strtod also accepts (case insensitively) strings of the form
138  *      NaN(x), where x is a string of hexadecimal digits and spaces;
139  *      if there is only one string of hexadecimal digits, it is taken
140  *      for the 52 fraction bits of the resulting NaN; if there are two
141  *      or more strings of hex digits, the first is for the high 20 bits,
142  *      the second and subsequent for the low 32 bits, with intervening
143  *      white space ignored; but if this results in none of the 52
144  *      fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
145  *      and NAN_WORD1 are used instead.
146  * #define MULTIPLE_THREADS if the system offers preemptively scheduled
147  *      multiple threads.  In this case, you must provide (or suitably
148  *      #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
149  *      by FREE_DTOA_LOCK(n) for n = 0 or 1.  (The second lock, accessed
150  *      in pow5mult, ensures lazy evaluation of only one copy of high
151  *      powers of 5; omitting this lock would introduce a small
152  *      probability of wasting memory, but would otherwise be harmless.)
153  *      You must also invoke freedtoa(s) to free the value s returned by
154  *      dtoa.  You may do so whether or not MULTIPLE_THREADS is #defined.
155
156  *      When MULTIPLE_THREADS is #defined, this source file provides
157  *              void set_max_dtoa_threads(unsigned int n);
158  *      and expects
159  *              unsigned int dtoa_get_threadno(void);
160  *      to be available (possibly provided by
161  *              #define dtoa_get_threadno omp_get_thread_num
162  *      if OpenMP is in use or by
163  *              #define dtoa_get_threadno pthread_self
164  *      if Pthreads is in use), to return the current thread number.
165  *      If set_max_dtoa_threads(n) was called and the current thread
166  *      number is k with k < n, then calls on ACQUIRE_DTOA_LOCK(...) and
167  *      FREE_DTOA_LOCK(...) are avoided; instead each thread with thread
168  *      number < n has a separate copy of relevant data structures.
169  *      After set_max_dtoa_threads(n), a call set_max_dtoa_threads(m)
170  *      with m <= n has has no effect, but a call with m > n is honored.
171  *      Such a call invokes REALLOC (assumed to be "realloc" if REALLOC
172  *      is not #defined) to extend the size of the relevant array.
173
174  * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
175  *      avoids underflows on inputs whose result does not underflow.
176  *      If you #define NO_IEEE_Scale on a machine that uses IEEE-format
177  *      floating-point numbers and flushes underflows to zero rather
178  *      than implementing gradual underflow, then you must also #define
179  *      Sudden_Underflow.
180  * #define USE_LOCALE to use the current locale's decimal_point value.
181  * #define SET_INEXACT if IEEE arithmetic is being used and extra
182  *      computation should be done to set the inexact flag when the
183  *      result is inexact and avoid setting inexact when the result
184  *      is exact.  In this case, dtoa.c must be compiled in
185  *      an environment, perhaps provided by #include "dtoa.c" in a
186  *      suitable wrapper, that defines two functions,
187  *              int get_inexact(void);
188  *              void clear_inexact(void);
189  *      such that get_inexact() returns a nonzero value if the
190  *      inexact bit is already set, and clear_inexact() sets the
191  *      inexact bit to 0.  When SET_INEXACT is #defined, strtod
192  *      also does extra computations to set the underflow and overflow
193  *      flags when appropriate (i.e., when the result is tiny and
194  *      inexact or when it is a numeric value rounded to +-infinity).
195  * #define NO_ERRNO if strtod should not assign errno = ERANGE when
196  *      the result overflows to +-Infinity or underflows to 0.
197  *      When errno should be assigned, under seemingly rare conditions
198  *      it may be necessary to define Set_errno(x) suitably, e.g., in
199  *      a local errno.h, such as
200  *              #include <errno.h>
201  *              #define Set_errno(x) _set_errno(x)
202  * #define NO_HEX_FP to omit recognition of hexadecimal floating-point
203  *      values by strtod.
204  * #define NO_STRTOD_BIGCOMP (on IEEE-arithmetic systems only for now)
205  *      to disable logic for "fast" testing of very long input strings
206  *      to strtod.  This testing proceeds by initially truncating the
207  *      input string, then if necessary comparing the whole string with
208  *      a decimal expansion to decide close cases. This logic is only
209  *      used for input more than STRTOD_DIGLIM digits long (default 40).
210  */
211
212 #ifndef Long
213 #define Long int
214 #endif
215 #ifndef ULong
216 typedef unsigned Long ULong;
217 #endif
218
219 #ifdef DEBUG
220 #include <assert.h>
221 #include "stdio.h"
222 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
223 #define Debug(x) x
224 int dtoa_stats[7]; /* strtod_{64,96,bigcomp},dtoa_{exact,64,96,bigcomp} */
225 #else
226 #define assert(x) /*nothing*/
227 #define Debug(x) /*nothing*/
228 #endif
229
230 #include "stdlib.h"
231 #include "string.h"
232
233 #ifdef USE_LOCALE
234 #include "locale.h"
235 #endif
236
237 #ifdef Honor_FLT_ROUNDS
238 #ifndef Trust_FLT_ROUNDS
239 #include <fenv.h>
240 #endif
241 #endif
242
243 #ifdef __cplusplus
244 extern "C" {
245 #endif
246 #ifdef MALLOC
247 extern void *MALLOC(size_t);
248 #else
249 #define MALLOC malloc
250 #endif
251
252 #ifdef REALLOC
253 extern void *REALLOC(void*,size_t);
254 #else
255 #define REALLOC realloc
256 #endif
257
258 #ifndef FREE
259 #define FREE free
260 #endif
261
262 #ifdef __cplusplus
263         }
264 #endif
265
266 #ifndef Omit_Private_Memory
267 #ifndef PRIVATE_MEM
268 #define PRIVATE_MEM 2304
269 #endif
270 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
271 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
272 #endif
273
274 #undef IEEE_Arith
275 #undef Avoid_Underflow
276 #ifdef IEEE_MC68k
277 #define IEEE_Arith
278 #endif
279 #ifdef IEEE_8087
280 #define IEEE_Arith
281 #endif
282
283 #ifdef IEEE_Arith
284 #ifndef NO_INFNAN_CHECK
285 #undef INFNAN_CHECK
286 #define INFNAN_CHECK
287 #endif
288 #else
289 #undef INFNAN_CHECK
290 #define NO_STRTOD_BIGCOMP
291 #endif
292
293 #include "errno.h"
294
295 #ifdef NO_ERRNO /*{*/
296 #undef Set_errno
297 #define Set_errno(x)
298 #else
299 #ifndef Set_errno
300 #define Set_errno(x) errno = x
301 #endif
302 #endif /*}*/
303
304 #ifdef Bad_float_h
305
306 #ifdef IEEE_Arith
307 #define DBL_DIG 15
308 #define DBL_MAX_10_EXP 308
309 #define DBL_MAX_EXP 1024
310 #define FLT_RADIX 2
311 #endif /*IEEE_Arith*/
312
313 #ifdef IBM
314 #define DBL_DIG 16
315 #define DBL_MAX_10_EXP 75
316 #define DBL_MAX_EXP 63
317 #define FLT_RADIX 16
318 #define DBL_MAX 7.2370055773322621e+75
319 #endif
320
321 #ifdef VAX
322 #define DBL_DIG 16
323 #define DBL_MAX_10_EXP 38
324 #define DBL_MAX_EXP 127
325 #define FLT_RADIX 2
326 #define DBL_MAX 1.7014118346046923e+38
327 #endif
328
329 #ifndef LONG_MAX
330 #define LONG_MAX 2147483647
331 #endif
332
333 #else /* ifndef Bad_float_h */
334 #include "float.h"
335 #endif /* Bad_float_h */
336
337 #ifndef __MATH_H__
338 #include "math.h"
339 #endif
340
341 #ifdef __cplusplus
342 extern "C" {
343 #endif
344
345 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
346 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
347 #endif
348
349 #undef USE_BF96
350
351 #ifdef NO_LONG_LONG /*{{*/
352 #undef ULLong
353 #ifdef Just_16
354 #undef Pack_32
355 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
356  * This makes some inner loops simpler and sometimes saves work
357  * during multiplications, but it often seems to make things slightly
358  * slower.  Hence the default is now to store 32 bits per Long.
359  */
360 #endif
361 #else   /*}{ long long available */
362 #ifndef Llong
363 #define Llong long long
364 #endif
365 #ifndef ULLong
366 #define ULLong unsigned Llong
367 #endif
368 #ifndef NO_BF96 /*{*/
369 #define USE_BF96
370
371 #ifdef SET_INEXACT
372 #define dtoa_divmax 27
373 #else
374 int dtoa_divmax = 2;    /* Permit experimenting: on some systems, 64-bit integer */
375                         /* division is slow enough that we may sometimes want to */
376                         /* avoid using it.   We assume (but do not check) that   */
377                         /* dtoa_divmax <= 27.*/
378 #endif
379
380 typedef struct BF96 {           /* Normalized 96-bit software floating point numbers */
381         unsigned int b0,b1,b2;  /* b0 = most significant, binary point just to its left */
382         int e;                  /* number represented = b * 2^e, with .5 <= b < 1 */
383         } BF96;
384
385  static BF96 pten[667] = {
386         { 0xeef453d6, 0x923bd65a, 0x113faa29, -1136 },
387         { 0x9558b466, 0x1b6565f8, 0x4ac7ca59, -1132 },
388         { 0xbaaee17f, 0xa23ebf76, 0x5d79bcf0, -1129 },
389         { 0xe95a99df, 0x8ace6f53, 0xf4d82c2c, -1126 },
390         { 0x91d8a02b, 0xb6c10594, 0x79071b9b, -1122 },
391         { 0xb64ec836, 0xa47146f9, 0x9748e282, -1119 },
392         { 0xe3e27a44, 0x4d8d98b7, 0xfd1b1b23, -1116 },
393         { 0x8e6d8c6a, 0xb0787f72, 0xfe30f0f5, -1112 },
394         { 0xb208ef85, 0x5c969f4f, 0xbdbd2d33, -1109 },
395         { 0xde8b2b66, 0xb3bc4723, 0xad2c7880, -1106 },
396         { 0x8b16fb20, 0x3055ac76, 0x4c3bcb50, -1102 },
397         { 0xaddcb9e8, 0x3c6b1793, 0xdf4abe24, -1099 },
398         { 0xd953e862, 0x4b85dd78, 0xd71d6dad, -1096 },
399         { 0x87d4713d, 0x6f33aa6b, 0x8672648c, -1092 },
400         { 0xa9c98d8c, 0xcb009506, 0x680efdaf, -1089 },
401         { 0xd43bf0ef, 0xfdc0ba48, 0x0212bd1b, -1086 },
402         { 0x84a57695, 0xfe98746d, 0x014bb630, -1082 },
403         { 0xa5ced43b, 0x7e3e9188, 0x419ea3bd, -1079 },
404         { 0xcf42894a, 0x5dce35ea, 0x52064cac, -1076 },
405         { 0x818995ce, 0x7aa0e1b2, 0x7343efeb, -1072 },
406         { 0xa1ebfb42, 0x19491a1f, 0x1014ebe6, -1069 },
407         { 0xca66fa12, 0x9f9b60a6, 0xd41a26e0, -1066 },
408         { 0xfd00b897, 0x478238d0, 0x8920b098, -1063 },
409         { 0x9e20735e, 0x8cb16382, 0x55b46e5f, -1059 },
410         { 0xc5a89036, 0x2fddbc62, 0xeb2189f7, -1056 },
411         { 0xf712b443, 0xbbd52b7b, 0xa5e9ec75, -1053 },
412         { 0x9a6bb0aa, 0x55653b2d, 0x47b233c9, -1049 },
413         { 0xc1069cd4, 0xeabe89f8, 0x999ec0bb, -1046 },
414         { 0xf148440a, 0x256e2c76, 0xc00670ea, -1043 },
415         { 0x96cd2a86, 0x5764dbca, 0x38040692, -1039 },
416         { 0xbc807527, 0xed3e12bc, 0xc6050837, -1036 },
417         { 0xeba09271, 0xe88d976b, 0xf7864a44, -1033 },
418         { 0x93445b87, 0x31587ea3, 0x7ab3ee6a, -1029 },
419         { 0xb8157268, 0xfdae9e4c, 0x5960ea05, -1026 },
420         { 0xe61acf03, 0x3d1a45df, 0x6fb92487, -1023 },
421         { 0x8fd0c162, 0x06306bab, 0xa5d3b6d4, -1019 },
422         { 0xb3c4f1ba, 0x87bc8696, 0x8f48a489, -1016 },
423         { 0xe0b62e29, 0x29aba83c, 0x331acdab, -1013 },
424         { 0x8c71dcd9, 0xba0b4925, 0x9ff0c08b, -1009 },
425         { 0xaf8e5410, 0x288e1b6f, 0x07ecf0ae, -1006 },
426         { 0xdb71e914, 0x32b1a24a, 0xc9e82cd9, -1003 },
427         { 0x892731ac, 0x9faf056e, 0xbe311c08,  -999 },
428         { 0xab70fe17, 0xc79ac6ca, 0x6dbd630a,  -996 },
429         { 0xd64d3d9d, 0xb981787d, 0x092cbbcc,  -993 },
430         { 0x85f04682, 0x93f0eb4e, 0x25bbf560,  -989 },
431         { 0xa76c5823, 0x38ed2621, 0xaf2af2b8,  -986 },
432         { 0xd1476e2c, 0x07286faa, 0x1af5af66,  -983 },
433         { 0x82cca4db, 0x847945ca, 0x50d98d9f,  -979 },
434         { 0xa37fce12, 0x6597973c, 0xe50ff107,  -976 },
435         { 0xcc5fc196, 0xfefd7d0c, 0x1e53ed49,  -973 },
436         { 0xff77b1fc, 0xbebcdc4f, 0x25e8e89c,  -970 },
437         { 0x9faacf3d, 0xf73609b1, 0x77b19161,  -966 },
438         { 0xc795830d, 0x75038c1d, 0xd59df5b9,  -963 },
439         { 0xf97ae3d0, 0xd2446f25, 0x4b057328,  -960 },
440         { 0x9becce62, 0x836ac577, 0x4ee367f9,  -956 },
441         { 0xc2e801fb, 0x244576d5, 0x229c41f7,  -953 },
442         { 0xf3a20279, 0xed56d48a, 0x6b435275,  -950 },
443         { 0x9845418c, 0x345644d6, 0x830a1389,  -946 },
444         { 0xbe5691ef, 0x416bd60c, 0x23cc986b,  -943 },
445         { 0xedec366b, 0x11c6cb8f, 0x2cbfbe86,  -940 },
446         { 0x94b3a202, 0xeb1c3f39, 0x7bf7d714,  -936 },
447         { 0xb9e08a83, 0xa5e34f07, 0xdaf5ccd9,  -933 },
448         { 0xe858ad24, 0x8f5c22c9, 0xd1b3400f,  -930 },
449         { 0x91376c36, 0xd99995be, 0x23100809,  -926 },
450         { 0xb5854744, 0x8ffffb2d, 0xabd40a0c,  -923 },
451         { 0xe2e69915, 0xb3fff9f9, 0x16c90c8f,  -920 },
452         { 0x8dd01fad, 0x907ffc3b, 0xae3da7d9,  -916 },
453         { 0xb1442798, 0xf49ffb4a, 0x99cd11cf,  -913 },
454         { 0xdd95317f, 0x31c7fa1d, 0x40405643,  -910 },
455         { 0x8a7d3eef, 0x7f1cfc52, 0x482835ea,  -906 },
456         { 0xad1c8eab, 0x5ee43b66, 0xda324365,  -903 },
457         { 0xd863b256, 0x369d4a40, 0x90bed43e,  -900 },
458         { 0x873e4f75, 0xe2224e68, 0x5a7744a6,  -896 },
459         { 0xa90de353, 0x5aaae202, 0x711515d0,  -893 },
460         { 0xd3515c28, 0x31559a83, 0x0d5a5b44,  -890 },
461         { 0x8412d999, 0x1ed58091, 0xe858790a,  -886 },
462         { 0xa5178fff, 0x668ae0b6, 0x626e974d,  -883 },
463         { 0xce5d73ff, 0x402d98e3, 0xfb0a3d21,  -880 },
464         { 0x80fa687f, 0x881c7f8e, 0x7ce66634,  -876 },
465         { 0xa139029f, 0x6a239f72, 0x1c1fffc1,  -873 },
466         { 0xc9874347, 0x44ac874e, 0xa327ffb2,  -870 },
467         { 0xfbe91419, 0x15d7a922, 0x4bf1ff9f,  -867 },
468         { 0x9d71ac8f, 0xada6c9b5, 0x6f773fc3,  -863 },
469         { 0xc4ce17b3, 0x99107c22, 0xcb550fb4,  -860 },
470         { 0xf6019da0, 0x7f549b2b, 0x7e2a53a1,  -857 },
471         { 0x99c10284, 0x4f94e0fb, 0x2eda7444,  -853 },
472         { 0xc0314325, 0x637a1939, 0xfa911155,  -850 },
473         { 0xf03d93ee, 0xbc589f88, 0x793555ab,  -847 },
474         { 0x96267c75, 0x35b763b5, 0x4bc1558b,  -843 },
475         { 0xbbb01b92, 0x83253ca2, 0x9eb1aaed,  -840 },
476         { 0xea9c2277, 0x23ee8bcb, 0x465e15a9,  -837 },
477         { 0x92a1958a, 0x7675175f, 0x0bfacd89,  -833 },
478         { 0xb749faed, 0x14125d36, 0xcef980ec,  -830 },
479         { 0xe51c79a8, 0x5916f484, 0x82b7e127,  -827 },
480         { 0x8f31cc09, 0x37ae58d2, 0xd1b2ecb8,  -823 },
481         { 0xb2fe3f0b, 0x8599ef07, 0x861fa7e6,  -820 },
482         { 0xdfbdcece, 0x67006ac9, 0x67a791e0,  -817 },
483         { 0x8bd6a141, 0x006042bd, 0xe0c8bb2c,  -813 },
484         { 0xaecc4991, 0x4078536d, 0x58fae9f7,  -810 },
485         { 0xda7f5bf5, 0x90966848, 0xaf39a475,  -807 },
486         { 0x888f9979, 0x7a5e012d, 0x6d8406c9,  -803 },
487         { 0xaab37fd7, 0xd8f58178, 0xc8e5087b,  -800 },
488         { 0xd5605fcd, 0xcf32e1d6, 0xfb1e4a9a,  -797 },
489         { 0x855c3be0, 0xa17fcd26, 0x5cf2eea0,  -793 },
490         { 0xa6b34ad8, 0xc9dfc06f, 0xf42faa48,  -790 },
491         { 0xd0601d8e, 0xfc57b08b, 0xf13b94da,  -787 },
492         { 0x823c1279, 0x5db6ce57, 0x76c53d08,  -783 },
493         { 0xa2cb1717, 0xb52481ed, 0x54768c4b,  -780 },
494         { 0xcb7ddcdd, 0xa26da268, 0xa9942f5d,  -777 },
495         { 0xfe5d5415, 0x0b090b02, 0xd3f93b35,  -774 },
496         { 0x9efa548d, 0x26e5a6e1, 0xc47bc501,  -770 },
497         { 0xc6b8e9b0, 0x709f109a, 0x359ab641,  -767 },
498         { 0xf867241c, 0x8cc6d4c0, 0xc30163d2,  -764 },
499         { 0x9b407691, 0xd7fc44f8, 0x79e0de63,  -760 },
500         { 0xc2109436, 0x4dfb5636, 0x985915fc,  -757 },
501         { 0xf294b943, 0xe17a2bc4, 0x3e6f5b7b,  -754 },
502         { 0x979cf3ca, 0x6cec5b5a, 0xa705992c,  -750 },
503         { 0xbd8430bd, 0x08277231, 0x50c6ff78,  -747 },
504         { 0xece53cec, 0x4a314ebd, 0xa4f8bf56,  -744 },
505         { 0x940f4613, 0xae5ed136, 0x871b7795,  -740 },
506         { 0xb9131798, 0x99f68584, 0x28e2557b,  -737 },
507         { 0xe757dd7e, 0xc07426e5, 0x331aeada,  -734 },
508         { 0x9096ea6f, 0x3848984f, 0x3ff0d2c8,  -730 },
509         { 0xb4bca50b, 0x065abe63, 0x0fed077a,  -727 },
510         { 0xe1ebce4d, 0xc7f16dfb, 0xd3e84959,  -724 },
511         { 0x8d3360f0, 0x9cf6e4bd, 0x64712dd7,  -720 },
512         { 0xb080392c, 0xc4349dec, 0xbd8d794d,  -717 },
513         { 0xdca04777, 0xf541c567, 0xecf0d7a0,  -714 },
514         { 0x89e42caa, 0xf9491b60, 0xf41686c4,  -710 },
515         { 0xac5d37d5, 0xb79b6239, 0x311c2875,  -707 },
516         { 0xd77485cb, 0x25823ac7, 0x7d633293,  -704 },
517         { 0x86a8d39e, 0xf77164bc, 0xae5dff9c,  -700 },
518         { 0xa8530886, 0xb54dbdeb, 0xd9f57f83,  -697 },
519         { 0xd267caa8, 0x62a12d66, 0xd072df63,  -694 },
520         { 0x8380dea9, 0x3da4bc60, 0x4247cb9e,  -690 },
521         { 0xa4611653, 0x8d0deb78, 0x52d9be85,  -687 },
522         { 0xcd795be8, 0x70516656, 0x67902e27,  -684 },
523         { 0x806bd971, 0x4632dff6, 0x00ba1cd8,  -680 },
524         { 0xa086cfcd, 0x97bf97f3, 0x80e8a40e,  -677 },
525         { 0xc8a883c0, 0xfdaf7df0, 0x6122cd12,  -674 },
526         { 0xfad2a4b1, 0x3d1b5d6c, 0x796b8057,  -671 },
527         { 0x9cc3a6ee, 0xc6311a63, 0xcbe33036,  -667 },
528         { 0xc3f490aa, 0x77bd60fc, 0xbedbfc44,  -664 },
529         { 0xf4f1b4d5, 0x15acb93b, 0xee92fb55,  -661 },
530         { 0x99171105, 0x2d8bf3c5, 0x751bdd15,  -657 },
531         { 0xbf5cd546, 0x78eef0b6, 0xd262d45a,  -654 },
532         { 0xef340a98, 0x172aace4, 0x86fb8971,  -651 },
533         { 0x9580869f, 0x0e7aac0e, 0xd45d35e6,  -647 },
534         { 0xbae0a846, 0xd2195712, 0x89748360,  -644 },
535         { 0xe998d258, 0x869facd7, 0x2bd1a438,  -641 },
536         { 0x91ff8377, 0x5423cc06, 0x7b6306a3,  -637 },
537         { 0xb67f6455, 0x292cbf08, 0x1a3bc84c,  -634 },
538         { 0xe41f3d6a, 0x7377eeca, 0x20caba5f,  -631 },
539         { 0x8e938662, 0x882af53e, 0x547eb47b,  -627 },
540         { 0xb23867fb, 0x2a35b28d, 0xe99e619a,  -624 },
541         { 0xdec681f9, 0xf4c31f31, 0x6405fa00,  -621 },
542         { 0x8b3c113c, 0x38f9f37e, 0xde83bc40,  -617 },
543         { 0xae0b158b, 0x4738705e, 0x9624ab50,  -614 },
544         { 0xd98ddaee, 0x19068c76, 0x3badd624,  -611 },
545         { 0x87f8a8d4, 0xcfa417c9, 0xe54ca5d7,  -607 },
546         { 0xa9f6d30a, 0x038d1dbc, 0x5e9fcf4c,  -604 },
547         { 0xd47487cc, 0x8470652b, 0x7647c320,  -601 },
548         { 0x84c8d4df, 0xd2c63f3b, 0x29ecd9f4,  -597 },
549         { 0xa5fb0a17, 0xc777cf09, 0xf4681071,  -594 },
550         { 0xcf79cc9d, 0xb955c2cc, 0x7182148d,  -591 },
551         { 0x81ac1fe2, 0x93d599bf, 0xc6f14cd8,  -587 },
552         { 0xa21727db, 0x38cb002f, 0xb8ada00e,  -584 },
553         { 0xca9cf1d2, 0x06fdc03b, 0xa6d90811,  -581 },
554         { 0xfd442e46, 0x88bd304a, 0x908f4a16,  -578 },
555         { 0x9e4a9cec, 0x15763e2e, 0x9a598e4e,  -574 },
556         { 0xc5dd4427, 0x1ad3cdba, 0x40eff1e1,  -571 },
557         { 0xf7549530, 0xe188c128, 0xd12bee59,  -568 },
558         { 0x9a94dd3e, 0x8cf578b9, 0x82bb74f8,  -564 },
559         { 0xc13a148e, 0x3032d6e7, 0xe36a5236,  -561 },
560         { 0xf18899b1, 0xbc3f8ca1, 0xdc44e6c3,  -558 },
561         { 0x96f5600f, 0x15a7b7e5, 0x29ab103a,  -554 },
562         { 0xbcb2b812, 0xdb11a5de, 0x7415d448,  -551 },
563         { 0xebdf6617, 0x91d60f56, 0x111b495b,  -548 },
564         { 0x936b9fce, 0xbb25c995, 0xcab10dd9,  -544 },
565         { 0xb84687c2, 0x69ef3bfb, 0x3d5d514f,  -541 },
566         { 0xe65829b3, 0x046b0afa, 0x0cb4a5a3,  -538 },
567         { 0x8ff71a0f, 0xe2c2e6dc, 0x47f0e785,  -534 },
568         { 0xb3f4e093, 0xdb73a093, 0x59ed2167,  -531 },
569         { 0xe0f218b8, 0xd25088b8, 0x306869c1,  -528 },
570         { 0x8c974f73, 0x83725573, 0x1e414218,  -524 },
571         { 0xafbd2350, 0x644eeacf, 0xe5d1929e,  -521 },
572         { 0xdbac6c24, 0x7d62a583, 0xdf45f746,  -518 },
573         { 0x894bc396, 0xce5da772, 0x6b8bba8c,  -514 },
574         { 0xab9eb47c, 0x81f5114f, 0x066ea92f,  -511 },
575         { 0xd686619b, 0xa27255a2, 0xc80a537b,  -508 },
576         { 0x8613fd01, 0x45877585, 0xbd06742c,  -504 },
577         { 0xa798fc41, 0x96e952e7, 0x2c481138,  -501 },
578         { 0xd17f3b51, 0xfca3a7a0, 0xf75a1586,  -498 },
579         { 0x82ef8513, 0x3de648c4, 0x9a984d73,  -494 },
580         { 0xa3ab6658, 0x0d5fdaf5, 0xc13e60d0,  -491 },
581         { 0xcc963fee, 0x10b7d1b3, 0x318df905,  -488 },
582         { 0xffbbcfe9, 0x94e5c61f, 0xfdf17746,  -485 },
583         { 0x9fd561f1, 0xfd0f9bd3, 0xfeb6ea8b,  -481 },
584         { 0xc7caba6e, 0x7c5382c8, 0xfe64a52e,  -478 },
585         { 0xf9bd690a, 0x1b68637b, 0x3dfdce7a,  -475 },
586         { 0x9c1661a6, 0x51213e2d, 0x06bea10c,  -471 },
587         { 0xc31bfa0f, 0xe5698db8, 0x486e494f,  -468 },
588         { 0xf3e2f893, 0xdec3f126, 0x5a89dba3,  -465 },
589         { 0x986ddb5c, 0x6b3a76b7, 0xf8962946,  -461 },
590         { 0xbe895233, 0x86091465, 0xf6bbb397,  -458 },
591         { 0xee2ba6c0, 0x678b597f, 0x746aa07d,  -455 },
592         { 0x94db4838, 0x40b717ef, 0xa8c2a44e,  -451 },
593         { 0xba121a46, 0x50e4ddeb, 0x92f34d62,  -448 },
594         { 0xe896a0d7, 0xe51e1566, 0x77b020ba,  -445 },
595         { 0x915e2486, 0xef32cd60, 0x0ace1474,  -441 },
596         { 0xb5b5ada8, 0xaaff80b8, 0x0d819992,  -438 },
597         { 0xe3231912, 0xd5bf60e6, 0x10e1fff6,  -435 },
598         { 0x8df5efab, 0xc5979c8f, 0xca8d3ffa,  -431 },
599         { 0xb1736b96, 0xb6fd83b3, 0xbd308ff8,  -428 },
600         { 0xddd0467c, 0x64bce4a0, 0xac7cb3f6,  -425 },
601         { 0x8aa22c0d, 0xbef60ee4, 0x6bcdf07a,  -421 },
602         { 0xad4ab711, 0x2eb3929d, 0x86c16c98,  -418 },
603         { 0xd89d64d5, 0x7a607744, 0xe871c7bf,  -415 },
604         { 0x87625f05, 0x6c7c4a8b, 0x11471cd7,  -411 },
605         { 0xa93af6c6, 0xc79b5d2d, 0xd598e40d,  -408 },
606         { 0xd389b478, 0x79823479, 0x4aff1d10,  -405 },
607         { 0x843610cb, 0x4bf160cb, 0xcedf722a,  -401 },
608         { 0xa54394fe, 0x1eedb8fe, 0xc2974eb4,  -398 },
609         { 0xce947a3d, 0xa6a9273e, 0x733d2262,  -395 },
610         { 0x811ccc66, 0x8829b887, 0x0806357d,  -391 },
611         { 0xa163ff80, 0x2a3426a8, 0xca07c2dc,  -388 },
612         { 0xc9bcff60, 0x34c13052, 0xfc89b393,  -385 },
613         { 0xfc2c3f38, 0x41f17c67, 0xbbac2078,  -382 },
614         { 0x9d9ba783, 0x2936edc0, 0xd54b944b,  -378 },
615         { 0xc5029163, 0xf384a931, 0x0a9e795e,  -375 },
616         { 0xf64335bc, 0xf065d37d, 0x4d4617b5,  -372 },
617         { 0x99ea0196, 0x163fa42e, 0x504bced1,  -368 },
618         { 0xc06481fb, 0x9bcf8d39, 0xe45ec286,  -365 },
619         { 0xf07da27a, 0x82c37088, 0x5d767327,  -362 },
620         { 0x964e858c, 0x91ba2655, 0x3a6a07f8,  -358 },
621         { 0xbbe226ef, 0xb628afea, 0x890489f7,  -355 },
622         { 0xeadab0ab, 0xa3b2dbe5, 0x2b45ac74,  -352 },
623         { 0x92c8ae6b, 0x464fc96f, 0x3b0b8bc9,  -348 },
624         { 0xb77ada06, 0x17e3bbcb, 0x09ce6ebb,  -345 },
625         { 0xe5599087, 0x9ddcaabd, 0xcc420a6a,  -342 },
626         { 0x8f57fa54, 0xc2a9eab6, 0x9fa94682,  -338 },
627         { 0xb32df8e9, 0xf3546564, 0x47939822,  -335 },
628         { 0xdff97724, 0x70297ebd, 0x59787e2b,  -332 },
629         { 0x8bfbea76, 0xc619ef36, 0x57eb4edb,  -328 },
630         { 0xaefae514, 0x77a06b03, 0xede62292,  -325 },
631         { 0xdab99e59, 0x958885c4, 0xe95fab36,  -322 },
632         { 0x88b402f7, 0xfd75539b, 0x11dbcb02,  -318 },
633         { 0xaae103b5, 0xfcd2a881, 0xd652bdc2,  -315 },
634         { 0xd59944a3, 0x7c0752a2, 0x4be76d33,  -312 },
635         { 0x857fcae6, 0x2d8493a5, 0x6f70a440,  -308 },
636         { 0xa6dfbd9f, 0xb8e5b88e, 0xcb4ccd50,  -305 },
637         { 0xd097ad07, 0xa71f26b2, 0x7e2000a4,  -302 },
638         { 0x825ecc24, 0xc873782f, 0x8ed40066,  -298 },
639         { 0xa2f67f2d, 0xfa90563b, 0x72890080,  -295 },
640         { 0xcbb41ef9, 0x79346bca, 0x4f2b40a0,  -292 },
641         { 0xfea126b7, 0xd78186bc, 0xe2f610c8,  -289 },
642         { 0x9f24b832, 0xe6b0f436, 0x0dd9ca7d,  -285 },
643         { 0xc6ede63f, 0xa05d3143, 0x91503d1c,  -282 },
644         { 0xf8a95fcf, 0x88747d94, 0x75a44c63,  -279 },
645         { 0x9b69dbe1, 0xb548ce7c, 0xc986afbe,  -275 },
646         { 0xc24452da, 0x229b021b, 0xfbe85bad,  -272 },
647         { 0xf2d56790, 0xab41c2a2, 0xfae27299,  -269 },
648         { 0x97c560ba, 0x6b0919a5, 0xdccd879f,  -265 },
649         { 0xbdb6b8e9, 0x05cb600f, 0x5400e987,  -262 },
650         { 0xed246723, 0x473e3813, 0x290123e9,  -259 },
651         { 0x9436c076, 0x0c86e30b, 0xf9a0b672,  -255 },
652         { 0xb9447093, 0x8fa89bce, 0xf808e40e,  -252 },
653         { 0xe7958cb8, 0x7392c2c2, 0xb60b1d12,  -249 },
654         { 0x90bd77f3, 0x483bb9b9, 0xb1c6f22b,  -245 },
655         { 0xb4ecd5f0, 0x1a4aa828, 0x1e38aeb6,  -242 },
656         { 0xe2280b6c, 0x20dd5232, 0x25c6da63,  -239 },
657         { 0x8d590723, 0x948a535f, 0x579c487e,  -235 },
658         { 0xb0af48ec, 0x79ace837, 0x2d835a9d,  -232 },
659         { 0xdcdb1b27, 0x98182244, 0xf8e43145,  -229 },
660         { 0x8a08f0f8, 0xbf0f156b, 0x1b8e9ecb,  -225 },
661         { 0xac8b2d36, 0xeed2dac5, 0xe272467e,  -222 },
662         { 0xd7adf884, 0xaa879177, 0x5b0ed81d,  -219 },
663         { 0x86ccbb52, 0xea94baea, 0x98e94712,  -215 },
664         { 0xa87fea27, 0xa539e9a5, 0x3f2398d7,  -212 },
665         { 0xd29fe4b1, 0x8e88640e, 0x8eec7f0d,  -209 },
666         { 0x83a3eeee, 0xf9153e89, 0x1953cf68,  -205 },
667         { 0xa48ceaaa, 0xb75a8e2b, 0x5fa8c342,  -202 },
668         { 0xcdb02555, 0x653131b6, 0x3792f412,  -199 },
669         { 0x808e1755, 0x5f3ebf11, 0xe2bbd88b,  -195 },
670         { 0xa0b19d2a, 0xb70e6ed6, 0x5b6aceae,  -192 },
671         { 0xc8de0475, 0x64d20a8b, 0xf245825a,  -189 },
672         { 0xfb158592, 0xbe068d2e, 0xeed6e2f0,  -186 },
673         { 0x9ced737b, 0xb6c4183d, 0x55464dd6,  -182 },
674         { 0xc428d05a, 0xa4751e4c, 0xaa97e14c,  -179 },
675         { 0xf5330471, 0x4d9265df, 0xd53dd99f,  -176 },
676         { 0x993fe2c6, 0xd07b7fab, 0xe546a803,  -172 },
677         { 0xbf8fdb78, 0x849a5f96, 0xde985204,  -169 },
678         { 0xef73d256, 0xa5c0f77c, 0x963e6685,  -166 },
679         { 0x95a86376, 0x27989aad, 0xdde70013,  -162 },
680         { 0xbb127c53, 0xb17ec159, 0x5560c018,  -159 },
681         { 0xe9d71b68, 0x9dde71af, 0xaab8f01e,  -156 },
682         { 0x92267121, 0x62ab070d, 0xcab39613,  -152 },
683         { 0xb6b00d69, 0xbb55c8d1, 0x3d607b97,  -149 },
684         { 0xe45c10c4, 0x2a2b3b05, 0x8cb89a7d,  -146 },
685         { 0x8eb98a7a, 0x9a5b04e3, 0x77f3608e,  -142 },
686         { 0xb267ed19, 0x40f1c61c, 0x55f038b2,  -139 },
687         { 0xdf01e85f, 0x912e37a3, 0x6b6c46de,  -136 },
688         { 0x8b61313b, 0xbabce2c6, 0x2323ac4b,  -132 },
689         { 0xae397d8a, 0xa96c1b77, 0xabec975e,  -129 },
690         { 0xd9c7dced, 0x53c72255, 0x96e7bd35,  -126 },
691         { 0x881cea14, 0x545c7575, 0x7e50d641,  -122 },
692         { 0xaa242499, 0x697392d2, 0xdde50bd1,  -119 },
693         { 0xd4ad2dbf, 0xc3d07787, 0x955e4ec6,  -116 },
694         { 0x84ec3c97, 0xda624ab4, 0xbd5af13b,  -112 },
695         { 0xa6274bbd, 0xd0fadd61, 0xecb1ad8a,  -109 },
696         { 0xcfb11ead, 0x453994ba, 0x67de18ed,  -106 },
697         { 0x81ceb32c, 0x4b43fcf4, 0x80eacf94,  -102 },
698         { 0xa2425ff7, 0x5e14fc31, 0xa1258379,   -99 },
699         { 0xcad2f7f5, 0x359a3b3e, 0x096ee458,   -96 },
700         { 0xfd87b5f2, 0x8300ca0d, 0x8bca9d6e,   -93 },
701         { 0x9e74d1b7, 0x91e07e48, 0x775ea264,   -89 },
702         { 0xc6120625, 0x76589dda, 0x95364afe,   -86 },
703         { 0xf79687ae, 0xd3eec551, 0x3a83ddbd,   -83 },
704         { 0x9abe14cd, 0x44753b52, 0xc4926a96,   -79 },
705         { 0xc16d9a00, 0x95928a27, 0x75b7053c,   -76 },
706         { 0xf1c90080, 0xbaf72cb1, 0x5324c68b,   -73 },
707         { 0x971da050, 0x74da7bee, 0xd3f6fc16,   -69 },
708         { 0xbce50864, 0x92111aea, 0x88f4bb1c,   -66 },
709         { 0xec1e4a7d, 0xb69561a5, 0x2b31e9e3,   -63 },
710         { 0x9392ee8e, 0x921d5d07, 0x3aff322e,   -59 },
711         { 0xb877aa32, 0x36a4b449, 0x09befeb9,   -56 },
712         { 0xe69594be, 0xc44de15b, 0x4c2ebe68,   -53 },
713         { 0x901d7cf7, 0x3ab0acd9, 0x0f9d3701,   -49 },
714         { 0xb424dc35, 0x095cd80f, 0x538484c1,   -46 },
715         { 0xe12e1342, 0x4bb40e13, 0x2865a5f2,   -43 },
716         { 0x8cbccc09, 0x6f5088cb, 0xf93f87b7,   -39 },
717         { 0xafebff0b, 0xcb24aafe, 0xf78f69a5,   -36 },
718         { 0xdbe6fece, 0xbdedd5be, 0xb573440e,   -33 },
719         { 0x89705f41, 0x36b4a597, 0x31680a88,   -29 },
720         { 0xabcc7711, 0x8461cefc, 0xfdc20d2b,   -26 },
721         { 0xd6bf94d5, 0xe57a42bc, 0x3d329076,   -23 },
722         { 0x8637bd05, 0xaf6c69b5, 0xa63f9a49,   -19 },
723         { 0xa7c5ac47, 0x1b478423, 0x0fcf80dc,   -16 },
724         { 0xd1b71758, 0xe219652b, 0xd3c36113,   -13 },
725         { 0x83126e97, 0x8d4fdf3b, 0x645a1cac,    -9 },
726         { 0xa3d70a3d, 0x70a3d70a, 0x3d70a3d7,    -6 },
727         { 0xcccccccc, 0xcccccccc, 0xcccccccc,    -3 },
728         { 0x80000000, 0x00000000, 0x00000000,    1 },
729         { 0xa0000000, 0x00000000, 0x00000000,    4 },
730         { 0xc8000000, 0x00000000, 0x00000000,    7 },
731         { 0xfa000000, 0x00000000, 0x00000000,   10 },
732         { 0x9c400000, 0x00000000, 0x00000000,   14 },
733         { 0xc3500000, 0x00000000, 0x00000000,   17 },
734         { 0xf4240000, 0x00000000, 0x00000000,   20 },
735         { 0x98968000, 0x00000000, 0x00000000,   24 },
736         { 0xbebc2000, 0x00000000, 0x00000000,   27 },
737         { 0xee6b2800, 0x00000000, 0x00000000,   30 },
738         { 0x9502f900, 0x00000000, 0x00000000,   34 },
739         { 0xba43b740, 0x00000000, 0x00000000,   37 },
740         { 0xe8d4a510, 0x00000000, 0x00000000,   40 },
741         { 0x9184e72a, 0x00000000, 0x00000000,   44 },
742         { 0xb5e620f4, 0x80000000, 0x00000000,   47 },
743         { 0xe35fa931, 0xa0000000, 0x00000000,   50 },
744         { 0x8e1bc9bf, 0x04000000, 0x00000000,   54 },
745         { 0xb1a2bc2e, 0xc5000000, 0x00000000,   57 },
746         { 0xde0b6b3a, 0x76400000, 0x00000000,   60 },
747         { 0x8ac72304, 0x89e80000, 0x00000000,   64 },
748         { 0xad78ebc5, 0xac620000, 0x00000000,   67 },
749         { 0xd8d726b7, 0x177a8000, 0x00000000,   70 },
750         { 0x87867832, 0x6eac9000, 0x00000000,   74 },
751         { 0xa968163f, 0x0a57b400, 0x00000000,   77 },
752         { 0xd3c21bce, 0xcceda100, 0x00000000,   80 },
753         { 0x84595161, 0x401484a0, 0x00000000,   84 },
754         { 0xa56fa5b9, 0x9019a5c8, 0x00000000,   87 },
755         { 0xcecb8f27, 0xf4200f3a, 0x00000000,   90 },
756         { 0x813f3978, 0xf8940984, 0x40000000,   94 },
757         { 0xa18f07d7, 0x36b90be5, 0x50000000,   97 },
758         { 0xc9f2c9cd, 0x04674ede, 0xa4000000,  100 },
759         { 0xfc6f7c40, 0x45812296, 0x4d000000,  103 },
760         { 0x9dc5ada8, 0x2b70b59d, 0xf0200000,  107 },
761         { 0xc5371912, 0x364ce305, 0x6c280000,  110 },
762         { 0xf684df56, 0xc3e01bc6, 0xc7320000,  113 },
763         { 0x9a130b96, 0x3a6c115c, 0x3c7f4000,  117 },
764         { 0xc097ce7b, 0xc90715b3, 0x4b9f1000,  120 },
765         { 0xf0bdc21a, 0xbb48db20, 0x1e86d400,  123 },
766         { 0x96769950, 0xb50d88f4, 0x13144480,  127 },
767         { 0xbc143fa4, 0xe250eb31, 0x17d955a0,  130 },
768         { 0xeb194f8e, 0x1ae525fd, 0x5dcfab08,  133 },
769         { 0x92efd1b8, 0xd0cf37be, 0x5aa1cae5,  137 },
770         { 0xb7abc627, 0x050305ad, 0xf14a3d9e,  140 },
771         { 0xe596b7b0, 0xc643c719, 0x6d9ccd05,  143 },
772         { 0x8f7e32ce, 0x7bea5c6f, 0xe4820023,  147 },
773         { 0xb35dbf82, 0x1ae4f38b, 0xdda2802c,  150 },
774         { 0xe0352f62, 0xa19e306e, 0xd50b2037,  153 },
775         { 0x8c213d9d, 0xa502de45, 0x4526f422,  157 },
776         { 0xaf298d05, 0x0e4395d6, 0x9670b12b,  160 },
777         { 0xdaf3f046, 0x51d47b4c, 0x3c0cdd76,  163 },
778         { 0x88d8762b, 0xf324cd0f, 0xa5880a69,  167 },
779         { 0xab0e93b6, 0xefee0053, 0x8eea0d04,  170 },
780         { 0xd5d238a4, 0xabe98068, 0x72a49045,  173 },
781         { 0x85a36366, 0xeb71f041, 0x47a6da2b,  177 },
782         { 0xa70c3c40, 0xa64e6c51, 0x999090b6,  180 },
783         { 0xd0cf4b50, 0xcfe20765, 0xfff4b4e3,  183 },
784         { 0x82818f12, 0x81ed449f, 0xbff8f10e,  187 },
785         { 0xa321f2d7, 0x226895c7, 0xaff72d52,  190 },
786         { 0xcbea6f8c, 0xeb02bb39, 0x9bf4f8a6,  193 },
787         { 0xfee50b70, 0x25c36a08, 0x02f236d0,  196 },
788         { 0x9f4f2726, 0x179a2245, 0x01d76242,  200 },
789         { 0xc722f0ef, 0x9d80aad6, 0x424d3ad2,  203 },
790         { 0xf8ebad2b, 0x84e0d58b, 0xd2e08987,  206 },
791         { 0x9b934c3b, 0x330c8577, 0x63cc55f4,  210 },
792         { 0xc2781f49, 0xffcfa6d5, 0x3cbf6b71,  213 },
793         { 0xf316271c, 0x7fc3908a, 0x8bef464e,  216 },
794         { 0x97edd871, 0xcfda3a56, 0x97758bf0,  220 },
795         { 0xbde94e8e, 0x43d0c8ec, 0x3d52eeed,  223 },
796         { 0xed63a231, 0xd4c4fb27, 0x4ca7aaa8,  226 },
797         { 0x945e455f, 0x24fb1cf8, 0x8fe8caa9,  230 },
798         { 0xb975d6b6, 0xee39e436, 0xb3e2fd53,  233 },
799         { 0xe7d34c64, 0xa9c85d44, 0x60dbbca8,  236 },
800         { 0x90e40fbe, 0xea1d3a4a, 0xbc8955e9,  240 },
801         { 0xb51d13ae, 0xa4a488dd, 0x6babab63,  243 },
802         { 0xe264589a, 0x4dcdab14, 0xc696963c,  246 },
803         { 0x8d7eb760, 0x70a08aec, 0xfc1e1de5,  250 },
804         { 0xb0de6538, 0x8cc8ada8, 0x3b25a55f,  253 },
805         { 0xdd15fe86, 0xaffad912, 0x49ef0eb7,  256 },
806         { 0x8a2dbf14, 0x2dfcc7ab, 0x6e356932,  260 },
807         { 0xacb92ed9, 0x397bf996, 0x49c2c37f,  263 },
808         { 0xd7e77a8f, 0x87daf7fb, 0xdc33745e,  266 },
809         { 0x86f0ac99, 0xb4e8dafd, 0x69a028bb,  270 },
810         { 0xa8acd7c0, 0x222311bc, 0xc40832ea,  273 },
811         { 0xd2d80db0, 0x2aabd62b, 0xf50a3fa4,  276 },
812         { 0x83c7088e, 0x1aab65db, 0x792667c6,  280 },
813         { 0xa4b8cab1, 0xa1563f52, 0x577001b8,  283 },
814         { 0xcde6fd5e, 0x09abcf26, 0xed4c0226,  286 },
815         { 0x80b05e5a, 0xc60b6178, 0x544f8158,  290 },
816         { 0xa0dc75f1, 0x778e39d6, 0x696361ae,  293 },
817         { 0xc913936d, 0xd571c84c, 0x03bc3a19,  296 },
818         { 0xfb587849, 0x4ace3a5f, 0x04ab48a0,  299 },
819         { 0x9d174b2d, 0xcec0e47b, 0x62eb0d64,  303 },
820         { 0xc45d1df9, 0x42711d9a, 0x3ba5d0bd,  306 },
821         { 0xf5746577, 0x930d6500, 0xca8f44ec,  309 },
822         { 0x9968bf6a, 0xbbe85f20, 0x7e998b13,  313 },
823         { 0xbfc2ef45, 0x6ae276e8, 0x9e3fedd8,  316 },
824         { 0xefb3ab16, 0xc59b14a2, 0xc5cfe94e,  319 },
825         { 0x95d04aee, 0x3b80ece5, 0xbba1f1d1,  323 },
826         { 0xbb445da9, 0xca61281f, 0x2a8a6e45,  326 },
827         { 0xea157514, 0x3cf97226, 0xf52d09d7,  329 },
828         { 0x924d692c, 0xa61be758, 0x593c2626,  333 },
829         { 0xb6e0c377, 0xcfa2e12e, 0x6f8b2fb0,  336 },
830         { 0xe498f455, 0xc38b997a, 0x0b6dfb9c,  339 },
831         { 0x8edf98b5, 0x9a373fec, 0x4724bd41,  343 },
832         { 0xb2977ee3, 0x00c50fe7, 0x58edec91,  346 },
833         { 0xdf3d5e9b, 0xc0f653e1, 0x2f2967b6,  349 },
834         { 0x8b865b21, 0x5899f46c, 0xbd79e0d2,  353 },
835         { 0xae67f1e9, 0xaec07187, 0xecd85906,  356 },
836         { 0xda01ee64, 0x1a708de9, 0xe80e6f48,  359 },
837         { 0x884134fe, 0x908658b2, 0x3109058d,  363 },
838         { 0xaa51823e, 0x34a7eede, 0xbd4b46f0,  366 },
839         { 0xd4e5e2cd, 0xc1d1ea96, 0x6c9e18ac,  369 },
840         { 0x850fadc0, 0x9923329e, 0x03e2cf6b,  373 },
841         { 0xa6539930, 0xbf6bff45, 0x84db8346,  376 },
842         { 0xcfe87f7c, 0xef46ff16, 0xe6126418,  379 },
843         { 0x81f14fae, 0x158c5f6e, 0x4fcb7e8f,  383 },
844         { 0xa26da399, 0x9aef7749, 0xe3be5e33,  386 },
845         { 0xcb090c80, 0x01ab551c, 0x5cadf5bf,  389 },
846         { 0xfdcb4fa0, 0x02162a63, 0x73d9732f,  392 },
847         { 0x9e9f11c4, 0x014dda7e, 0x2867e7fd,  396 },
848         { 0xc646d635, 0x01a1511d, 0xb281e1fd,  399 },
849         { 0xf7d88bc2, 0x4209a565, 0x1f225a7c,  402 },
850         { 0x9ae75759, 0x6946075f, 0x3375788d,  406 },
851         { 0xc1a12d2f, 0xc3978937, 0x0052d6b1,  409 },
852         { 0xf209787b, 0xb47d6b84, 0xc0678c5d,  412 },
853         { 0x9745eb4d, 0x50ce6332, 0xf840b7ba,  416 },
854         { 0xbd176620, 0xa501fbff, 0xb650e5a9,  419 },
855         { 0xec5d3fa8, 0xce427aff, 0xa3e51f13,  422 },
856         { 0x93ba47c9, 0x80e98cdf, 0xc66f336c,  426 },
857         { 0xb8a8d9bb, 0xe123f017, 0xb80b0047,  429 },
858         { 0xe6d3102a, 0xd96cec1d, 0xa60dc059,  432 },
859         { 0x9043ea1a, 0xc7e41392, 0x87c89837,  436 },
860         { 0xb454e4a1, 0x79dd1877, 0x29babe45,  439 },
861         { 0xe16a1dc9, 0xd8545e94, 0xf4296dd6,  442 },
862         { 0x8ce2529e, 0x2734bb1d, 0x1899e4a6,  446 },
863         { 0xb01ae745, 0xb101e9e4, 0x5ec05dcf,  449 },
864         { 0xdc21a117, 0x1d42645d, 0x76707543,  452 },
865         { 0x899504ae, 0x72497eba, 0x6a06494a,  456 },
866         { 0xabfa45da, 0x0edbde69, 0x0487db9d,  459 },
867         { 0xd6f8d750, 0x9292d603, 0x45a9d284,  462 },
868         { 0x865b8692, 0x5b9bc5c2, 0x0b8a2392,  466 },
869         { 0xa7f26836, 0xf282b732, 0x8e6cac77,  469 },
870         { 0xd1ef0244, 0xaf2364ff, 0x3207d795,  472 },
871         { 0x8335616a, 0xed761f1f, 0x7f44e6bd,  476 },
872         { 0xa402b9c5, 0xa8d3a6e7, 0x5f16206c,  479 },
873         { 0xcd036837, 0x130890a1, 0x36dba887,  482 },
874         { 0x80222122, 0x6be55a64, 0xc2494954,  486 },
875         { 0xa02aa96b, 0x06deb0fd, 0xf2db9baa,  489 },
876         { 0xc83553c5, 0xc8965d3d, 0x6f928294,  492 },
877         { 0xfa42a8b7, 0x3abbf48c, 0xcb772339,  495 },
878         { 0x9c69a972, 0x84b578d7, 0xff2a7604,  499 },
879         { 0xc38413cf, 0x25e2d70d, 0xfef51385,  502 },
880         { 0xf46518c2, 0xef5b8cd1, 0x7eb25866,  505 },
881         { 0x98bf2f79, 0xd5993802, 0xef2f773f,  509 },
882         { 0xbeeefb58, 0x4aff8603, 0xaafb550f,  512 },
883         { 0xeeaaba2e, 0x5dbf6784, 0x95ba2a53,  515 },
884         { 0x952ab45c, 0xfa97a0b2, 0xdd945a74,  519 },
885         { 0xba756174, 0x393d88df, 0x94f97111,  522 },
886         { 0xe912b9d1, 0x478ceb17, 0x7a37cd56,  525 },
887         { 0x91abb422, 0xccb812ee, 0xac62e055,  529 },
888         { 0xb616a12b, 0x7fe617aa, 0x577b986b,  532 },
889         { 0xe39c4976, 0x5fdf9d94, 0xed5a7e85,  535 },
890         { 0x8e41ade9, 0xfbebc27d, 0x14588f13,  539 },
891         { 0xb1d21964, 0x7ae6b31c, 0x596eb2d8,  542 },
892         { 0xde469fbd, 0x99a05fe3, 0x6fca5f8e,  545 },
893         { 0x8aec23d6, 0x80043bee, 0x25de7bb9,  549 },
894         { 0xada72ccc, 0x20054ae9, 0xaf561aa7,  552 },
895         { 0xd910f7ff, 0x28069da4, 0x1b2ba151,  555 },
896         { 0x87aa9aff, 0x79042286, 0x90fb44d2,  559 },
897         { 0xa99541bf, 0x57452b28, 0x353a1607,  562 },
898         { 0xd3fa922f, 0x2d1675f2, 0x42889b89,  565 },
899         { 0x847c9b5d, 0x7c2e09b7, 0x69956135,  569 },
900         { 0xa59bc234, 0xdb398c25, 0x43fab983,  572 },
901         { 0xcf02b2c2, 0x1207ef2e, 0x94f967e4,  575 },
902         { 0x8161afb9, 0x4b44f57d, 0x1d1be0ee,  579 },
903         { 0xa1ba1ba7, 0x9e1632dc, 0x6462d92a,  582 },
904         { 0xca28a291, 0x859bbf93, 0x7d7b8f75,  585 },
905         { 0xfcb2cb35, 0xe702af78, 0x5cda7352,  588 },
906         { 0x9defbf01, 0xb061adab, 0x3a088813,  592 },
907         { 0xc56baec2, 0x1c7a1916, 0x088aaa18,  595 },
908         { 0xf6c69a72, 0xa3989f5b, 0x8aad549e,  598 },
909         { 0x9a3c2087, 0xa63f6399, 0x36ac54e2,  602 },
910         { 0xc0cb28a9, 0x8fcf3c7f, 0x84576a1b,  605 },
911         { 0xf0fdf2d3, 0xf3c30b9f, 0x656d44a2,  608 },
912         { 0x969eb7c4, 0x7859e743, 0x9f644ae5,  612 },
913         { 0xbc4665b5, 0x96706114, 0x873d5d9f,  615 },
914         { 0xeb57ff22, 0xfc0c7959, 0xa90cb506,  618 },
915         { 0x9316ff75, 0xdd87cbd8, 0x09a7f124,  622 },
916         { 0xb7dcbf53, 0x54e9bece, 0x0c11ed6d,  625 },
917         { 0xe5d3ef28, 0x2a242e81, 0x8f1668c8,  628 },
918         { 0x8fa47579, 0x1a569d10, 0xf96e017d,  632 },
919         { 0xb38d92d7, 0x60ec4455, 0x37c981dc,  635 },
920         { 0xe070f78d, 0x3927556a, 0x85bbe253,  638 },
921         { 0x8c469ab8, 0x43b89562, 0x93956d74,  642 },
922         { 0xaf584166, 0x54a6babb, 0x387ac8d1,  645 },
923         { 0xdb2e51bf, 0xe9d0696a, 0x06997b05,  648 },
924         { 0x88fcf317, 0xf22241e2, 0x441fece3,  652 },
925         { 0xab3c2fdd, 0xeeaad25a, 0xd527e81c,  655 },
926         { 0xd60b3bd5, 0x6a5586f1, 0x8a71e223,  658 },
927         { 0x85c70565, 0x62757456, 0xf6872d56,  662 },
928         { 0xa738c6be, 0xbb12d16c, 0xb428f8ac,  665 },
929         { 0xd106f86e, 0x69d785c7, 0xe13336d7,  668 },
930         { 0x82a45b45, 0x0226b39c, 0xecc00246,  672 },
931         { 0xa34d7216, 0x42b06084, 0x27f002d7,  675 },
932         { 0xcc20ce9b, 0xd35c78a5, 0x31ec038d,  678 },
933         { 0xff290242, 0xc83396ce, 0x7e670471,  681 },
934         { 0x9f79a169, 0xbd203e41, 0x0f0062c6,  685 },
935         { 0xc75809c4, 0x2c684dd1, 0x52c07b78,  688 },
936         { 0xf92e0c35, 0x37826145, 0xa7709a56,  691 },
937         { 0x9bbcc7a1, 0x42b17ccb, 0x88a66076,  695 },
938         { 0xc2abf989, 0x935ddbfe, 0x6acff893,  698 },
939         { 0xf356f7eb, 0xf83552fe, 0x0583f6b8,  701 },
940         { 0x98165af3, 0x7b2153de, 0xc3727a33,  705 },
941         { 0xbe1bf1b0, 0x59e9a8d6, 0x744f18c0,  708 },
942         { 0xeda2ee1c, 0x7064130c, 0x1162def0,  711 },
943         { 0x9485d4d1, 0xc63e8be7, 0x8addcb56,  715 },
944         { 0xb9a74a06, 0x37ce2ee1, 0x6d953e2b,  718 },
945         { 0xe8111c87, 0xc5c1ba99, 0xc8fa8db6,  721 },
946         { 0x910ab1d4, 0xdb9914a0, 0x1d9c9892,  725 },
947         { 0xb54d5e4a, 0x127f59c8, 0x2503beb6,  728 },
948         { 0xe2a0b5dc, 0x971f303a, 0x2e44ae64,  731 },
949         { 0x8da471a9, 0xde737e24, 0x5ceaecfe,  735 },
950         { 0xb10d8e14, 0x56105dad, 0x7425a83e,  738 },
951         { 0xdd50f199, 0x6b947518, 0xd12f124e,  741 },
952         { 0x8a5296ff, 0xe33cc92f, 0x82bd6b70,  745 },
953         { 0xace73cbf, 0xdc0bfb7b, 0x636cc64d,  748 },
954         { 0xd8210bef, 0xd30efa5a, 0x3c47f7e0,  751 },
955         { 0x8714a775, 0xe3e95c78, 0x65acfaec,  755 },
956         { 0xa8d9d153, 0x5ce3b396, 0x7f1839a7,  758 },
957         { 0xd31045a8, 0x341ca07c, 0x1ede4811,  761 },
958         { 0x83ea2b89, 0x2091e44d, 0x934aed0a,  765 },
959         { 0xa4e4b66b, 0x68b65d60, 0xf81da84d,  768 },
960         { 0xce1de406, 0x42e3f4b9, 0x36251260,  771 },
961         { 0x80d2ae83, 0xe9ce78f3, 0xc1d72b7c,  775 },
962         { 0xa1075a24, 0xe4421730, 0xb24cf65b,  778 },
963         { 0xc94930ae, 0x1d529cfc, 0xdee033f2,  781 },
964         { 0xfb9b7cd9, 0xa4a7443c, 0x169840ef,  784 },
965         { 0x9d412e08, 0x06e88aa5, 0x8e1f2895,  788 },
966         { 0xc491798a, 0x08a2ad4e, 0xf1a6f2ba,  791 },
967         { 0xf5b5d7ec, 0x8acb58a2, 0xae10af69,  794 },
968         { 0x9991a6f3, 0xd6bf1765, 0xacca6da1,  798 },
969         { 0xbff610b0, 0xcc6edd3f, 0x17fd090a,  801 },
970         { 0xeff394dc, 0xff8a948e, 0xddfc4b4c,  804 },
971         { 0x95f83d0a, 0x1fb69cd9, 0x4abdaf10,  808 },
972         { 0xbb764c4c, 0xa7a4440f, 0x9d6d1ad4,  811 },
973         { 0xea53df5f, 0xd18d5513, 0x84c86189,  814 },
974         { 0x92746b9b, 0xe2f8552c, 0x32fd3cf5,  818 },
975         { 0xb7118682, 0xdbb66a77, 0x3fbc8c33,  821 },
976         { 0xe4d5e823, 0x92a40515, 0x0fabaf3f,  824 },
977         { 0x8f05b116, 0x3ba6832d, 0x29cb4d87,  828 },
978         { 0xb2c71d5b, 0xca9023f8, 0x743e20e9,  831 },
979         { 0xdf78e4b2, 0xbd342cf6, 0x914da924,  834 },
980         { 0x8bab8eef, 0xb6409c1a, 0x1ad089b6,  838 },
981         { 0xae9672ab, 0xa3d0c320, 0xa184ac24,  841 },
982         { 0xda3c0f56, 0x8cc4f3e8, 0xc9e5d72d,  844 },
983         { 0x88658996, 0x17fb1871, 0x7e2fa67c,  848 },
984         { 0xaa7eebfb, 0x9df9de8d, 0xddbb901b,  851 },
985         { 0xd51ea6fa, 0x85785631, 0x552a7422,  854 },
986         { 0x8533285c, 0x936b35de, 0xd53a8895,  858 },
987         { 0xa67ff273, 0xb8460356, 0x8a892aba,  861 },
988         { 0xd01fef10, 0xa657842c, 0x2d2b7569,  864 },
989         { 0x8213f56a, 0x67f6b29b, 0x9c3b2962,  868 },
990         { 0xa298f2c5, 0x01f45f42, 0x8349f3ba,  871 },
991         { 0xcb3f2f76, 0x42717713, 0x241c70a9,  874 },
992         { 0xfe0efb53, 0xd30dd4d7, 0xed238cd3,  877 },
993         { 0x9ec95d14, 0x63e8a506, 0xf4363804,  881 },
994         { 0xc67bb459, 0x7ce2ce48, 0xb143c605,  884 },
995         { 0xf81aa16f, 0xdc1b81da, 0xdd94b786,  887 },
996         { 0x9b10a4e5, 0xe9913128, 0xca7cf2b4,  891 },
997         { 0xc1d4ce1f, 0x63f57d72, 0xfd1c2f61,  894 },
998         { 0xf24a01a7, 0x3cf2dccf, 0xbc633b39,  897 },
999         { 0x976e4108, 0x8617ca01, 0xd5be0503,  901 },
1000         { 0xbd49d14a, 0xa79dbc82, 0x4b2d8644,  904 },
1001         { 0xec9c459d, 0x51852ba2, 0xddf8e7d6,  907 },
1002         { 0x93e1ab82, 0x52f33b45, 0xcabb90e5,  911 },
1003         { 0xb8da1662, 0xe7b00a17, 0x3d6a751f,  914 },
1004         { 0xe7109bfb, 0xa19c0c9d, 0x0cc51267,  917 },
1005         { 0x906a617d, 0x450187e2, 0x27fb2b80,  921 },
1006         { 0xb484f9dc, 0x9641e9da, 0xb1f9f660,  924 },
1007         { 0xe1a63853, 0xbbd26451, 0x5e7873f8,  927 },
1008         { 0x8d07e334, 0x55637eb2, 0xdb0b487b,  931 },
1009         { 0xb049dc01, 0x6abc5e5f, 0x91ce1a9a,  934 },
1010         { 0xdc5c5301, 0xc56b75f7, 0x7641a140,  937 },
1011         { 0x89b9b3e1, 0x1b6329ba, 0xa9e904c8,  941 },
1012         { 0xac2820d9, 0x623bf429, 0x546345fa,  944 },
1013         { 0xd732290f, 0xbacaf133, 0xa97c1779,  947 },
1014         { 0x867f59a9, 0xd4bed6c0, 0x49ed8eab,  951 },
1015         { 0xa81f3014, 0x49ee8c70, 0x5c68f256,  954 },
1016         { 0xd226fc19, 0x5c6a2f8c, 0x73832eec,  957 },
1017         { 0x83585d8f, 0xd9c25db7, 0xc831fd53,  961 },
1018         { 0xa42e74f3, 0xd032f525, 0xba3e7ca8,  964 },
1019         { 0xcd3a1230, 0xc43fb26f, 0x28ce1bd2,  967 },
1020         { 0x80444b5e, 0x7aa7cf85, 0x7980d163,  971 },
1021         { 0xa0555e36, 0x1951c366, 0xd7e105bc,  974 },
1022         { 0xc86ab5c3, 0x9fa63440, 0x8dd9472b,  977 },
1023         { 0xfa856334, 0x878fc150, 0xb14f98f6,  980 },
1024         { 0x9c935e00, 0xd4b9d8d2, 0x6ed1bf9a,  984 },
1025         { 0xc3b83581, 0x09e84f07, 0x0a862f80,  987 },
1026         { 0xf4a642e1, 0x4c6262c8, 0xcd27bb61,  990 },
1027         { 0x98e7e9cc, 0xcfbd7dbd, 0x8038d51c,  994 },
1028         { 0xbf21e440, 0x03acdd2c, 0xe0470a63,  997 },
1029         { 0xeeea5d50, 0x04981478, 0x1858ccfc, 1000 },
1030         { 0x95527a52, 0x02df0ccb, 0x0f37801e, 1004 },
1031         { 0xbaa718e6, 0x8396cffd, 0xd3056025, 1007 },
1032         { 0xe950df20, 0x247c83fd, 0x47c6b82e, 1010 },
1033         { 0x91d28b74, 0x16cdd27e, 0x4cdc331d, 1014 },
1034         { 0xb6472e51, 0x1c81471d, 0xe0133fe4, 1017 },
1035         { 0xe3d8f9e5, 0x63a198e5, 0x58180fdd, 1020 },
1036         { 0x8e679c2f, 0x5e44ff8f, 0x570f09ea, 1024 },
1037         { 0xb201833b, 0x35d63f73, 0x2cd2cc65, 1027 },
1038         { 0xde81e40a, 0x034bcf4f, 0xf8077f7e, 1030 },
1039         { 0x8b112e86, 0x420f6191, 0xfb04afaf, 1034 },
1040         { 0xadd57a27, 0xd29339f6, 0x79c5db9a, 1037 },
1041         { 0xd94ad8b1, 0xc7380874, 0x18375281, 1040 },
1042         { 0x87cec76f, 0x1c830548, 0x8f229391, 1044 },
1043         { 0xa9c2794a, 0xe3a3c69a, 0xb2eb3875, 1047 },
1044         { 0xd433179d, 0x9c8cb841, 0x5fa60692, 1050 },
1045         { 0x849feec2, 0x81d7f328, 0xdbc7c41b, 1054 },
1046         { 0xa5c7ea73, 0x224deff3, 0x12b9b522, 1057 },
1047         { 0xcf39e50f, 0xeae16bef, 0xd768226b, 1060 },
1048         { 0x81842f29, 0xf2cce375, 0xe6a11583, 1064 },
1049         { 0xa1e53af4, 0x6f801c53, 0x60495ae3, 1067 },
1050         { 0xca5e89b1, 0x8b602368, 0x385bb19c, 1070 },
1051         { 0xfcf62c1d, 0xee382c42, 0x46729e03, 1073 },
1052         { 0x9e19db92, 0xb4e31ba9, 0x6c07a2c2, 1077 }
1053         };
1054  static short int Lhint[2098] = {
1055            /*18,*/19,    19,    19,    19,    20,    20,    20,    21,    21,
1056            21,    22,    22,    22,    23,    23,    23,    23,    24,    24,
1057            24,    25,    25,    25,    26,    26,    26,    26,    27,    27,
1058            27,    28,    28,    28,    29,    29,    29,    29,    30,    30,
1059            30,    31,    31,    31,    32,    32,    32,    32,    33,    33,
1060            33,    34,    34,    34,    35,    35,    35,    35,    36,    36,
1061            36,    37,    37,    37,    38,    38,    38,    38,    39,    39,
1062            39,    40,    40,    40,    41,    41,    41,    41,    42,    42,
1063            42,    43,    43,    43,    44,    44,    44,    44,    45,    45,
1064            45,    46,    46,    46,    47,    47,    47,    47,    48,    48,
1065            48,    49,    49,    49,    50,    50,    50,    51,    51,    51,
1066            51,    52,    52,    52,    53,    53,    53,    54,    54,    54,
1067            54,    55,    55,    55,    56,    56,    56,    57,    57,    57,
1068            57,    58,    58,    58,    59,    59,    59,    60,    60,    60,
1069            60,    61,    61,    61,    62,    62,    62,    63,    63,    63,
1070            63,    64,    64,    64,    65,    65,    65,    66,    66,    66,
1071            66,    67,    67,    67,    68,    68,    68,    69,    69,    69,
1072            69,    70,    70,    70,    71,    71,    71,    72,    72,    72,
1073            72,    73,    73,    73,    74,    74,    74,    75,    75,    75,
1074            75,    76,    76,    76,    77,    77,    77,    78,    78,    78,
1075            78,    79,    79,    79,    80,    80,    80,    81,    81,    81,
1076            82,    82,    82,    82,    83,    83,    83,    84,    84,    84,
1077            85,    85,    85,    85,    86,    86,    86,    87,    87,    87,
1078            88,    88,    88,    88,    89,    89,    89,    90,    90,    90,
1079            91,    91,    91,    91,    92,    92,    92,    93,    93,    93,
1080            94,    94,    94,    94,    95,    95,    95,    96,    96,    96,
1081            97,    97,    97,    97,    98,    98,    98,    99,    99,    99,
1082           100,   100,   100,   100,   101,   101,   101,   102,   102,   102,
1083           103,   103,   103,   103,   104,   104,   104,   105,   105,   105,
1084           106,   106,   106,   106,   107,   107,   107,   108,   108,   108,
1085           109,   109,   109,   110,   110,   110,   110,   111,   111,   111,
1086           112,   112,   112,   113,   113,   113,   113,   114,   114,   114,
1087           115,   115,   115,   116,   116,   116,   116,   117,   117,   117,
1088           118,   118,   118,   119,   119,   119,   119,   120,   120,   120,
1089           121,   121,   121,   122,   122,   122,   122,   123,   123,   123,
1090           124,   124,   124,   125,   125,   125,   125,   126,   126,   126,
1091           127,   127,   127,   128,   128,   128,   128,   129,   129,   129,
1092           130,   130,   130,   131,   131,   131,   131,   132,   132,   132,
1093           133,   133,   133,   134,   134,   134,   134,   135,   135,   135,
1094           136,   136,   136,   137,   137,   137,   137,   138,   138,   138,
1095           139,   139,   139,   140,   140,   140,   141,   141,   141,   141,
1096           142,   142,   142,   143,   143,   143,   144,   144,   144,   144,
1097           145,   145,   145,   146,   146,   146,   147,   147,   147,   147,
1098           148,   148,   148,   149,   149,   149,   150,   150,   150,   150,
1099           151,   151,   151,   152,   152,   152,   153,   153,   153,   153,
1100           154,   154,   154,   155,   155,   155,   156,   156,   156,   156,
1101           157,   157,   157,   158,   158,   158,   159,   159,   159,   159,
1102           160,   160,   160,   161,   161,   161,   162,   162,   162,   162,
1103           163,   163,   163,   164,   164,   164,   165,   165,   165,   165,
1104           166,   166,   166,   167,   167,   167,   168,   168,   168,   169,
1105           169,   169,   169,   170,   170,   170,   171,   171,   171,   172,
1106           172,   172,   172,   173,   173,   173,   174,   174,   174,   175,
1107           175,   175,   175,   176,   176,   176,   177,   177,   177,   178,
1108           178,   178,   178,   179,   179,   179,   180,   180,   180,   181,
1109           181,   181,   181,   182,   182,   182,   183,   183,   183,   184,
1110           184,   184,   184,   185,   185,   185,   186,   186,   186,   187,
1111           187,   187,   187,   188,   188,   188,   189,   189,   189,   190,
1112           190,   190,   190,   191,   191,   191,   192,   192,   192,   193,
1113           193,   193,   193,   194,   194,   194,   195,   195,   195,   196,
1114           196,   196,   197,   197,   197,   197,   198,   198,   198,   199,
1115           199,   199,   200,   200,   200,   200,   201,   201,   201,   202,
1116           202,   202,   203,   203,   203,   203,   204,   204,   204,   205,
1117           205,   205,   206,   206,   206,   206,   207,   207,   207,   208,
1118           208,   208,   209,   209,   209,   209,   210,   210,   210,   211,
1119           211,   211,   212,   212,   212,   212,   213,   213,   213,   214,
1120           214,   214,   215,   215,   215,   215,   216,   216,   216,   217,
1121           217,   217,   218,   218,   218,   218,   219,   219,   219,   220,
1122           220,   220,   221,   221,   221,   221,   222,   222,   222,   223,
1123           223,   223,   224,   224,   224,   224,   225,   225,   225,   226,
1124           226,   226,   227,   227,   227,   228,   228,   228,   228,   229,
1125           229,   229,   230,   230,   230,   231,   231,   231,   231,   232,
1126           232,   232,   233,   233,   233,   234,   234,   234,   234,   235,
1127           235,   235,   236,   236,   236,   237,   237,   237,   237,   238,
1128           238,   238,   239,   239,   239,   240,   240,   240,   240,   241,
1129           241,   241,   242,   242,   242,   243,   243,   243,   243,   244,
1130           244,   244,   245,   245,   245,   246,   246,   246,   246,   247,
1131           247,   247,   248,   248,   248,   249,   249,   249,   249,   250,
1132           250,   250,   251,   251,   251,   252,   252,   252,   252,   253,
1133           253,   253,   254,   254,   254,   255,   255,   255,   256,   256,
1134           256,   256,   257,   257,   257,   258,   258,   258,   259,   259,
1135           259,   259,   260,   260,   260,   261,   261,   261,   262,   262,
1136           262,   262,   263,   263,   263,   264,   264,   264,   265,   265,
1137           265,   265,   266,   266,   266,   267,   267,   267,   268,   268,
1138           268,   268,   269,   269,   269,   270,   270,   270,   271,   271,
1139           271,   271,   272,   272,   272,   273,   273,   273,   274,   274,
1140           274,   274,   275,   275,   275,   276,   276,   276,   277,   277,
1141           277,   277,   278,   278,   278,   279,   279,   279,   280,   280,
1142           280,   280,   281,   281,   281,   282,   282,   282,   283,   283,
1143           283,   283,   284,   284,   284,   285,   285,   285,   286,   286,
1144           286,   287,   287,   287,   287,   288,   288,   288,   289,   289,
1145           289,   290,   290,   290,   290,   291,   291,   291,   292,   292,
1146           292,   293,   293,   293,   293,   294,   294,   294,   295,   295,
1147           295,   296,   296,   296,   296,   297,   297,   297,   298,   298,
1148           298,   299,   299,   299,   299,   300,   300,   300,   301,   301,
1149           301,   302,   302,   302,   302,   303,   303,   303,   304,   304,
1150           304,   305,   305,   305,   305,   306,   306,   306,   307,   307,
1151           307,   308,   308,   308,   308,   309,   309,   309,   310,   310,
1152           310,   311,   311,   311,   311,   312,   312,   312,   313,   313,
1153           313,   314,   314,   314,   315,   315,   315,   315,   316,   316,
1154           316,   317,   317,   317,   318,   318,   318,   318,   319,   319,
1155           319,   320,   320,   320,   321,   321,   321,   321,   322,   322,
1156           322,   323,   323,   323,   324,   324,   324,   324,   325,   325,
1157           325,   326,   326,   326,   327,   327,   327,   327,   328,   328,
1158           328,   329,   329,   329,   330,   330,   330,   330,   331,   331,
1159           331,   332,   332,   332,   333,   333,   333,   333,   334,   334,
1160           334,   335,   335,   335,   336,   336,   336,   336,   337,   337,
1161           337,   338,   338,   338,   339,   339,   339,   339,   340,   340,
1162           340,   341,   341,   341,   342,   342,   342,   342,   343,   343,
1163           343,   344,   344,   344,   345,   345,   345,   346,   346,   346,
1164           346,   347,   347,   347,   348,   348,   348,   349,   349,   349,
1165           349,   350,   350,   350,   351,   351,   351,   352,   352,   352,
1166           352,   353,   353,   353,   354,   354,   354,   355,   355,   355,
1167           355,   356,   356,   356,   357,   357,   357,   358,   358,   358,
1168           358,   359,   359,   359,   360,   360,   360,   361,   361,   361,
1169           361,   362,   362,   362,   363,   363,   363,   364,   364,   364,
1170           364,   365,   365,   365,   366,   366,   366,   367,   367,   367,
1171           367,   368,   368,   368,   369,   369,   369,   370,   370,   370,
1172           370,   371,   371,   371,   372,   372,   372,   373,   373,   373,
1173           374,   374,   374,   374,   375,   375,   375,   376,   376,   376,
1174           377,   377,   377,   377,   378,   378,   378,   379,   379,   379,
1175           380,   380,   380,   380,   381,   381,   381,   382,   382,   382,
1176           383,   383,   383,   383,   384,   384,   384,   385,   385,   385,
1177           386,   386,   386,   386,   387,   387,   387,   388,   388,   388,
1178           389,   389,   389,   389,   390,   390,   390,   391,   391,   391,
1179           392,   392,   392,   392,   393,   393,   393,   394,   394,   394,
1180           395,   395,   395,   395,   396,   396,   396,   397,   397,   397,
1181           398,   398,   398,   398,   399,   399,   399,   400,   400,   400,
1182           401,   401,   401,   402,   402,   402,   402,   403,   403,   403,
1183           404,   404,   404,   405,   405,   405,   405,   406,   406,   406,
1184           407,   407,   407,   408,   408,   408,   408,   409,   409,   409,
1185           410,   410,   410,   411,   411,   411,   411,   412,   412,   412,
1186           413,   413,   413,   414,   414,   414,   414,   415,   415,   415,
1187           416,   416,   416,   417,   417,   417,   417,   418,   418,   418,
1188           419,   419,   419,   420,   420,   420,   420,   421,   421,   421,
1189           422,   422,   422,   423,   423,   423,   423,   424,   424,   424,
1190           425,   425,   425,   426,   426,   426,   426,   427,   427,   427,
1191           428,   428,   428,   429,   429,   429,   429,   430,   430,   430,
1192           431,   431,   431,   432,   432,   432,   433,   433,   433,   433,
1193           434,   434,   434,   435,   435,   435,   436,   436,   436,   436,
1194           437,   437,   437,   438,   438,   438,   439,   439,   439,   439,
1195           440,   440,   440,   441,   441,   441,   442,   442,   442,   442,
1196           443,   443,   443,   444,   444,   444,   445,   445,   445,   445,
1197           446,   446,   446,   447,   447,   447,   448,   448,   448,   448,
1198           449,   449,   449,   450,   450,   450,   451,   451,   451,   451,
1199           452,   452,   452,   453,   453,   453,   454,   454,   454,   454,
1200           455,   455,   455,   456,   456,   456,   457,   457,   457,   457,
1201           458,   458,   458,   459,   459,   459,   460,   460,   460,   461,
1202           461,   461,   461,   462,   462,   462,   463,   463,   463,   464,
1203           464,   464,   464,   465,   465,   465,   466,   466,   466,   467,
1204           467,   467,   467,   468,   468,   468,   469,   469,   469,   470,
1205           470,   470,   470,   471,   471,   471,   472,   472,   472,   473,
1206           473,   473,   473,   474,   474,   474,   475,   475,   475,   476,
1207           476,   476,   476,   477,   477,   477,   478,   478,   478,   479,
1208           479,   479,   479,   480,   480,   480,   481,   481,   481,   482,
1209           482,   482,   482,   483,   483,   483,   484,   484,   484,   485,
1210           485,   485,   485,   486,   486,   486,   487,   487,   487,   488,
1211           488,   488,   488,   489,   489,   489,   490,   490,   490,   491,
1212           491,   491,   492,   492,   492,   492,   493,   493,   493,   494,
1213           494,   494,   495,   495,   495,   495,   496,   496,   496,   497,
1214           497,   497,   498,   498,   498,   498,   499,   499,   499,   500,
1215           500,   500,   501,   501,   501,   501,   502,   502,   502,   503,
1216           503,   503,   504,   504,   504,   504,   505,   505,   505,   506,
1217           506,   506,   507,   507,   507,   507,   508,   508,   508,   509,
1218           509,   509,   510,   510,   510,   510,   511,   511,   511,   512,
1219           512,   512,   513,   513,   513,   513,   514,   514,   514,   515,
1220           515,   515,   516,   516,   516,   516,   517,   517,   517,   518,
1221           518,   518,   519,   519,   519,   520,   520,   520,   520,   521,
1222           521,   521,   522,   522,   522,   523,   523,   523,   523,   524,
1223           524,   524,   525,   525,   525,   526,   526,   526,   526,   527,
1224           527,   527,   528,   528,   528,   529,   529,   529,   529,   530,
1225           530,   530,   531,   531,   531,   532,   532,   532,   532,   533,
1226           533,   533,   534,   534,   534,   535,   535,   535,   535,   536,
1227           536,   536,   537,   537,   537,   538,   538,   538,   538,   539,
1228           539,   539,   540,   540,   540,   541,   541,   541,   541,   542,
1229           542,   542,   543,   543,   543,   544,   544,   544,   544,   545,
1230           545,   545,   546,   546,   546,   547,   547,   547,   548,   548,
1231           548,   548,   549,   549,   549,   550,   550,   550,   551,   551,
1232           551,   551,   552,   552,   552,   553,   553,   553,   554,   554,
1233           554,   554,   555,   555,   555,   556,   556,   556,   557,   557,
1234           557,   557,   558,   558,   558,   559,   559,   559,   560,   560,
1235           560,   560,   561,   561,   561,   562,   562,   562,   563,   563,
1236           563,   563,   564,   564,   564,   565,   565,   565,   566,   566,
1237           566,   566,   567,   567,   567,   568,   568,   568,   569,   569,
1238           569,   569,   570,   570,   570,   571,   571,   571,   572,   572,
1239           572,   572,   573,   573,   573,   574,   574,   574,   575,   575,
1240           575,   575,   576,   576,   576,   577,   577,   577,   578,   578,
1241           578,   579,   579,   579,   579,   580,   580,   580,   581,   581,
1242           581,   582,   582,   582,   582,   583,   583,   583,   584,   584,
1243           584,   585,   585,   585,   585,   586,   586,   586,   587,   587,
1244           587,   588,   588,   588,   588,   589,   589,   589,   590,   590,
1245           590,   591,   591,   591,   591,   592,   592,   592,   593,   593,
1246           593,   594,   594,   594,   594,   595,   595,   595,   596,   596,
1247           596,   597,   597,   597,   597,   598,   598,   598,   599,   599,
1248           599,   600,   600,   600,   600,   601,   601,   601,   602,   602,
1249           602,   603,   603,   603,   603,   604,   604,   604,   605,   605,
1250           605,   606,   606,   606,   607,   607,   607,   607,   608,   608,
1251           608,   609,   609,   609,   610,   610,   610,   610,   611,   611,
1252           611,   612,   612,   612,   613,   613,   613,   613,   614,   614,
1253           614,   615,   615,   615,   616,   616,   616,   616,   617,   617,
1254           617,   618,   618,   618,   619,   619,   619,   619,   620,   620,
1255           620,   621,   621,   621,   622,   622,   622,   622,   623,   623,
1256           623,   624,   624,   624,   625,   625,   625,   625,   626,   626,
1257           626,   627,   627,   627,   628,   628,   628,   628,   629,   629,
1258           629,   630,   630,   630,   631,   631,   631,   631,   632,   632,
1259           632,   633,   633,   633,   634,   634,   634,   634,   635,   635,
1260           635,   636,   636,   636,   637,   637,   637,   638,   638,   638,
1261           638,   639,   639,   639,   640,   640,   640,   641,   641,   641,
1262           641,   642,   642,   642,   643,   643,   643,   644,   644,   644,
1263           644,   645,   645,   645,   646,   646,   646,   647,   647,   647,
1264           647,   648,   648,   648,   649,   649,   649,   650,   650 };
1265  static ULLong pfive[27] = {
1266                 5ll,
1267                 25ll,
1268                 125ll,
1269                 625ll,
1270                 3125ll,
1271                 15625ll,
1272                 78125ll,
1273                 390625ll,
1274                 1953125ll,
1275                 9765625ll,
1276                 48828125ll,
1277                 244140625ll,
1278                 1220703125ll,
1279                 6103515625ll,
1280                 30517578125ll,
1281                 152587890625ll,
1282                 762939453125ll,
1283                 3814697265625ll,
1284                 19073486328125ll,
1285                 95367431640625ll,
1286                 476837158203125ll,
1287                 2384185791015625ll,
1288                 11920928955078125ll,
1289                 59604644775390625ll,
1290                 298023223876953125ll,
1291                 1490116119384765625ll,
1292                 7450580596923828125ll
1293                 };
1294
1295  static int pfivebits[25] = {3, 5, 7, 10, 12, 14, 17, 19, 21, 24, 26, 28, 31,
1296                              33, 35, 38, 40, 42, 45, 47, 49, 52, 54, 56, 59};
1297 #endif /*}*/
1298 #endif /*}} NO_LONG_LONG */
1299
1300 typedef union { double d; ULong L[2];
1301 #ifdef USE_BF96
1302         ULLong LL;
1303 #endif
1304         } U;
1305
1306 #ifdef IEEE_8087
1307 #define word0(x) (x)->L[1]
1308 #define word1(x) (x)->L[0]
1309 #else
1310 #define word0(x) (x)->L[0]
1311 #define word1(x) (x)->L[1]
1312 #endif
1313 #define dval(x) (x)->d
1314 #define LLval(x) (x)->LL
1315
1316 #ifndef STRTOD_DIGLIM
1317 #define STRTOD_DIGLIM 40
1318 #endif
1319
1320 #ifdef DIGLIM_DEBUG
1321 extern int strtod_diglim;
1322 #else
1323 #define strtod_diglim STRTOD_DIGLIM
1324 #endif
1325
1326 /* The following definition of Storeinc is appropriate for MIPS processors.
1327  * An alternative that might be better on some machines is
1328  * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
1329  */
1330 #if defined(IEEE_8087) + defined(VAX)
1331 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
1332 ((unsigned short *)a)[0] = (unsigned short)c, a++)
1333 #else
1334 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
1335 ((unsigned short *)a)[1] = (unsigned short)c, a++)
1336 #endif
1337
1338 /* #define P DBL_MANT_DIG */
1339 /* Ten_pmax = floor(P*log(2)/log(5)) */
1340 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
1341 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
1342 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
1343
1344 #ifdef IEEE_Arith
1345 #define Exp_shift  20
1346 #define Exp_shift1 20
1347 #define Exp_msk1    0x100000
1348 #define Exp_msk11   0x100000
1349 #define Exp_mask  0x7ff00000
1350 #define P 53
1351 #define Nbits 53
1352 #define Bias 1023
1353 #define Emax 1023
1354 #define Emin (-1022)
1355 #define Exp_1  0x3ff00000
1356 #define Exp_11 0x3ff00000
1357 #define Ebits 11
1358 #define Frac_mask  0xfffff
1359 #define Frac_mask1 0xfffff
1360 #define Ten_pmax 22
1361 #define Bletch 0x10
1362 #define Bndry_mask  0xfffff
1363 #define Bndry_mask1 0xfffff
1364 #define LSB 1
1365 #define Sign_bit 0x80000000
1366 #define Log2P 1
1367 #define Tiny0 0
1368 #define Tiny1 1
1369 #define Quick_max 14
1370 #define Int_max 14
1371 #ifndef NO_IEEE_Scale
1372 #define Avoid_Underflow
1373 #ifdef Flush_Denorm     /* debugging option */
1374 #undef Sudden_Underflow
1375 #endif
1376 #endif
1377
1378 #ifndef Flt_Rounds
1379 #ifdef FLT_ROUNDS
1380 #define Flt_Rounds FLT_ROUNDS
1381 #else
1382 #define Flt_Rounds 1
1383 #endif
1384 #endif /*Flt_Rounds*/
1385
1386 #ifdef Honor_FLT_ROUNDS
1387 #undef Check_FLT_ROUNDS
1388 #define Check_FLT_ROUNDS
1389 #else
1390 #define Rounding Flt_Rounds
1391 #endif
1392
1393 #else /* ifndef IEEE_Arith */
1394 #undef Check_FLT_ROUNDS
1395 #undef Honor_FLT_ROUNDS
1396 #undef SET_INEXACT
1397 #undef  Sudden_Underflow
1398 #define Sudden_Underflow
1399 #ifdef IBM
1400 #undef Flt_Rounds
1401 #define Flt_Rounds 0
1402 #define Exp_shift  24
1403 #define Exp_shift1 24
1404 #define Exp_msk1   0x1000000
1405 #define Exp_msk11  0x1000000
1406 #define Exp_mask  0x7f000000
1407 #define P 14
1408 #define Nbits 56
1409 #define Bias 65
1410 #define Emax 248
1411 #define Emin (-260)
1412 #define Exp_1  0x41000000
1413 #define Exp_11 0x41000000
1414 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
1415 #define Frac_mask  0xffffff
1416 #define Frac_mask1 0xffffff
1417 #define Bletch 4
1418 #define Ten_pmax 22
1419 #define Bndry_mask  0xefffff
1420 #define Bndry_mask1 0xffffff
1421 #define LSB 1
1422 #define Sign_bit 0x80000000
1423 #define Log2P 4
1424 #define Tiny0 0x100000
1425 #define Tiny1 0
1426 #define Quick_max 14
1427 #define Int_max 15
1428 #else /* VAX */
1429 #undef Flt_Rounds
1430 #define Flt_Rounds 1
1431 #define Exp_shift  23
1432 #define Exp_shift1 7
1433 #define Exp_msk1    0x80
1434 #define Exp_msk11   0x800000
1435 #define Exp_mask  0x7f80
1436 #define P 56
1437 #define Nbits 56
1438 #define Bias 129
1439 #define Emax 126
1440 #define Emin (-129)
1441 #define Exp_1  0x40800000
1442 #define Exp_11 0x4080
1443 #define Ebits 8
1444 #define Frac_mask  0x7fffff
1445 #define Frac_mask1 0xffff007f
1446 #define Ten_pmax 24
1447 #define Bletch 2
1448 #define Bndry_mask  0xffff007f
1449 #define Bndry_mask1 0xffff007f
1450 #define LSB 0x10000
1451 #define Sign_bit 0x8000
1452 #define Log2P 1
1453 #define Tiny0 0x80
1454 #define Tiny1 0
1455 #define Quick_max 15
1456 #define Int_max 15
1457 #endif /* IBM, VAX */
1458 #endif /* IEEE_Arith */
1459
1460 #ifndef IEEE_Arith
1461 #define ROUND_BIASED
1462 #else
1463 #ifdef ROUND_BIASED_without_Round_Up
1464 #undef  ROUND_BIASED
1465 #define ROUND_BIASED
1466 #endif
1467 #endif
1468
1469 #ifdef RND_PRODQUOT
1470 #define rounded_product(a,b) a = rnd_prod(a, b)
1471 #define rounded_quotient(a,b) a = rnd_quot(a, b)
1472 extern double rnd_prod(double, double), rnd_quot(double, double);
1473 #else
1474 #define rounded_product(a,b) a *= b
1475 #define rounded_quotient(a,b) a /= b
1476 #endif
1477
1478 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
1479 #define Big1 0xffffffff
1480
1481 #ifndef Pack_32
1482 #define Pack_32
1483 #endif
1484
1485 typedef struct BCinfo BCinfo;
1486  struct
1487 BCinfo { int dp0, dp1, dplen, dsign, e0, inexact, nd, nd0, rounding, scale, uflchk; };
1488
1489 #define FFFFFFFF 0xffffffffUL
1490
1491 #ifdef MULTIPLE_THREADS
1492 #define MTa , PTI
1493 #define MTb , &TI
1494 #define MTd , ThInfo **PTI
1495 static unsigned int maxthreads = 0;
1496 #else
1497 #define MTa /*nothing*/
1498 #define MTb /*nothing*/
1499 #define MTd /*nothing*/
1500 #endif
1501
1502 #define Kmax 7
1503
1504 #ifdef __cplusplus
1505 extern "C" double strtod(const char *s00, char **se);
1506 extern "C" char *dtoa(double d, int mode, int ndigits,
1507                         int *decpt, int *sign, char **rve);
1508 #endif
1509
1510  struct
1511 Bigint {
1512         struct Bigint *next;
1513         int k, maxwds, sign, wds;
1514         ULong x[1];
1515         };
1516
1517  typedef struct Bigint Bigint;
1518  typedef struct
1519 ThInfo {
1520         Bigint *Freelist[Kmax+1];
1521         Bigint *P5s;
1522         } ThInfo;
1523
1524  static ThInfo TI0;
1525
1526 #ifdef MULTIPLE_THREADS
1527  static ThInfo *TI1;
1528  static int TI0_used;
1529
1530  void
1531 set_max_dtoa_threads(unsigned int n)
1532 {
1533         size_t L;
1534
1535         if (n > maxthreads) {
1536                 L = n*sizeof(ThInfo);
1537                 if (TI1) {
1538                         TI1 = (ThInfo*)REALLOC(TI1, L);
1539                         memset(TI1 + maxthreads, 0, (n-maxthreads)*sizeof(ThInfo));
1540                         }
1541                 else {
1542                         TI1 = (ThInfo*)MALLOC(L);
1543                         if (TI0_used) {
1544                                 memcpy(TI1, &TI0, sizeof(ThInfo));
1545                                 if (n > 1)
1546                                         memset(TI1 + 1, 0, L - sizeof(ThInfo));
1547                                 memset(&TI0, 0, sizeof(ThInfo));
1548                                 }
1549                         else
1550                                 memset(TI1, 0, L);
1551                         }
1552                 maxthreads = n;
1553                 }
1554         }
1555
1556  static ThInfo*
1557 get_TI(void)
1558 {
1559         unsigned int thno = dtoa_get_threadno();
1560         if (thno < maxthreads)
1561                 return TI1 + thno;
1562         if (thno == 0)
1563                 TI0_used = 1;
1564         return &TI0;
1565         }
1566 #define freelist TI->Freelist
1567 #define p5s TI->P5s
1568 #else
1569 #define freelist TI0.Freelist
1570 #define p5s TI0.P5s
1571 #endif
1572
1573  static Bigint *
1574 Balloc(int k MTd)
1575 {
1576         int x;
1577         Bigint *rv;
1578 #ifndef Omit_Private_Memory
1579         unsigned int len;
1580 #endif
1581 #ifdef MULTIPLE_THREADS
1582         ThInfo *TI;
1583
1584         if (!(TI = *PTI))
1585                 *PTI = TI = get_TI();
1586         if (TI == &TI0)
1587                 ACQUIRE_DTOA_LOCK(0);
1588 #endif
1589         /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
1590         /* but this case seems very unlikely. */
1591         if (k <= Kmax && (rv = freelist[k]))
1592                 freelist[k] = rv->next;
1593         else {
1594                 x = 1 << k;
1595 #ifdef Omit_Private_Memory
1596                 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
1597 #else
1598                 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
1599                         /sizeof(double);
1600                 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem
1601 #ifdef MULTIPLE_THREADS
1602                         && TI == TI1
1603 #endif
1604                         ) {
1605                         rv = (Bigint*)pmem_next;
1606                         pmem_next += len;
1607                         }
1608                 else
1609                         rv = (Bigint*)MALLOC(len*sizeof(double));
1610 #endif
1611                 rv->k = k;
1612                 rv->maxwds = x;
1613                 }
1614 #ifdef MULTIPLE_THREADS
1615         if (TI == &TI0)
1616                 FREE_DTOA_LOCK(0);
1617 #endif
1618         rv->sign = rv->wds = 0;
1619         return rv;
1620         }
1621
1622  static void
1623 Bfree(Bigint *v MTd)
1624 {
1625 #ifdef MULTIPLE_THREADS
1626         ThInfo *TI;
1627 #endif
1628         if (v) {
1629                 if (v->k > Kmax)
1630                         FREE((void*)v);
1631                 else {
1632 #ifdef MULTIPLE_THREADS
1633                         if (!(TI = *PTI))
1634                                 *PTI = TI = get_TI();
1635                         if (TI == &TI0)
1636                                 ACQUIRE_DTOA_LOCK(0);
1637 #endif
1638                         v->next = freelist[v->k];
1639                         freelist[v->k] = v;
1640 #ifdef MULTIPLE_THREADS
1641                         if (TI == &TI0)
1642                                 FREE_DTOA_LOCK(0);
1643 #endif
1644                         }
1645                 }
1646         }
1647
1648 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
1649 y->wds*sizeof(Long) + 2*sizeof(int))
1650
1651  static Bigint *
1652 multadd(Bigint *b, int m, int a MTd)    /* multiply by m and add a */
1653 {
1654         int i, wds;
1655 #ifdef ULLong
1656         ULong *x;
1657         ULLong carry, y;
1658 #else
1659         ULong carry, *x, y;
1660 #ifdef Pack_32
1661         ULong xi, z;
1662 #endif
1663 #endif
1664         Bigint *b1;
1665
1666         wds = b->wds;
1667         x = b->x;
1668         i = 0;
1669         carry = a;
1670         do {
1671 #ifdef ULLong
1672                 y = *x * (ULLong)m + carry;
1673                 carry = y >> 32;
1674                 *x++ = y & FFFFFFFF;
1675 #else
1676 #ifdef Pack_32
1677                 xi = *x;
1678                 y = (xi & 0xffff) * m + carry;
1679                 z = (xi >> 16) * m + (y >> 16);
1680                 carry = z >> 16;
1681                 *x++ = (z << 16) + (y & 0xffff);
1682 #else
1683                 y = *x * m + carry;
1684                 carry = y >> 16;
1685                 *x++ = y & 0xffff;
1686 #endif
1687 #endif
1688                 }
1689                 while(++i < wds);
1690         if (carry) {
1691                 if (wds >= b->maxwds) {
1692                         b1 = Balloc(b->k+1 MTa);
1693                         Bcopy(b1, b);
1694                         Bfree(b MTa);
1695                         b = b1;
1696                         }
1697                 b->x[wds++] = carry;
1698                 b->wds = wds;
1699                 }
1700         return b;
1701         }
1702
1703  static Bigint *
1704 s2b(const char *s, int nd0, int nd, ULong y9, int dplen MTd)
1705 {
1706         Bigint *b;
1707         int i, k;
1708         Long x, y;
1709
1710         x = (nd + 8) / 9;
1711         for(k = 0, y = 1; x > y; y <<= 1, k++) ;
1712 #ifdef Pack_32
1713         b = Balloc(k MTa);
1714         b->x[0] = y9;
1715         b->wds = 1;
1716 #else
1717         b = Balloc(k+1 MTa);
1718         b->x[0] = y9 & 0xffff;
1719         b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
1720 #endif
1721
1722         i = 9;
1723         if (9 < nd0) {
1724                 s += 9;
1725                 do b = multadd(b, 10, *s++ - '0' MTa);
1726                         while(++i < nd0);
1727                 s += dplen;
1728                 }
1729         else
1730                 s += dplen + 9;
1731         for(; i < nd; i++)
1732                 b = multadd(b, 10, *s++ - '0' MTa);
1733         return b;
1734         }
1735
1736  static int
1737 hi0bits(ULong x)
1738 {
1739         int k = 0;
1740
1741         if (!(x & 0xffff0000)) {
1742                 k = 16;
1743                 x <<= 16;
1744                 }
1745         if (!(x & 0xff000000)) {
1746                 k += 8;
1747                 x <<= 8;
1748                 }
1749         if (!(x & 0xf0000000)) {
1750                 k += 4;
1751                 x <<= 4;
1752                 }
1753         if (!(x & 0xc0000000)) {
1754                 k += 2;
1755                 x <<= 2;
1756                 }
1757         if (!(x & 0x80000000)) {
1758                 k++;
1759                 if (!(x & 0x40000000))
1760                         return 32;
1761                 }
1762         return k;
1763         }
1764
1765  static int
1766 lo0bits(ULong *y)
1767 {
1768         int k;
1769         ULong x = *y;
1770
1771         if (x & 7) {
1772                 if (x & 1)
1773                         return 0;
1774                 if (x & 2) {
1775                         *y = x >> 1;
1776                         return 1;
1777                         }
1778                 *y = x >> 2;
1779                 return 2;
1780                 }
1781         k = 0;
1782         if (!(x & 0xffff)) {
1783                 k = 16;
1784                 x >>= 16;
1785                 }
1786         if (!(x & 0xff)) {
1787                 k += 8;
1788                 x >>= 8;
1789                 }
1790         if (!(x & 0xf)) {
1791                 k += 4;
1792                 x >>= 4;
1793                 }
1794         if (!(x & 0x3)) {
1795                 k += 2;
1796                 x >>= 2;
1797                 }
1798         if (!(x & 1)) {
1799                 k++;
1800                 x >>= 1;
1801                 if (!x)
1802                         return 32;
1803                 }
1804         *y = x;
1805         return k;
1806         }
1807
1808  static Bigint *
1809 i2b(int i MTd)
1810 {
1811         Bigint *b;
1812
1813         b = Balloc(1 MTa);
1814         b->x[0] = i;
1815         b->wds = 1;
1816         return b;
1817         }
1818
1819  static Bigint *
1820 mult(Bigint *a, Bigint *b MTd)
1821 {
1822         Bigint *c;
1823         int k, wa, wb, wc;
1824         ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
1825         ULong y;
1826 #ifdef ULLong
1827         ULLong carry, z;
1828 #else
1829         ULong carry, z;
1830 #ifdef Pack_32
1831         ULong z2;
1832 #endif
1833 #endif
1834
1835         if (a->wds < b->wds) {
1836                 c = a;
1837                 a = b;
1838                 b = c;
1839                 }
1840         k = a->k;
1841         wa = a->wds;
1842         wb = b->wds;
1843         wc = wa + wb;
1844         if (wc > a->maxwds)
1845                 k++;
1846         c = Balloc(k MTa);
1847         for(x = c->x, xa = x + wc; x < xa; x++)
1848                 *x = 0;
1849         xa = a->x;
1850         xae = xa + wa;
1851         xb = b->x;
1852         xbe = xb + wb;
1853         xc0 = c->x;
1854 #ifdef ULLong
1855         for(; xb < xbe; xc0++) {
1856                 if ((y = *xb++)) {
1857                         x = xa;
1858                         xc = xc0;
1859                         carry = 0;
1860                         do {
1861                                 z = *x++ * (ULLong)y + *xc + carry;
1862                                 carry = z >> 32;
1863                                 *xc++ = z & FFFFFFFF;
1864                                 }
1865                                 while(x < xae);
1866                         *xc = carry;
1867                         }
1868                 }
1869 #else
1870 #ifdef Pack_32
1871         for(; xb < xbe; xb++, xc0++) {
1872                 if (y = *xb & 0xffff) {
1873                         x = xa;
1874                         xc = xc0;
1875                         carry = 0;
1876                         do {
1877                                 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
1878                                 carry = z >> 16;
1879                                 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
1880                                 carry = z2 >> 16;
1881                                 Storeinc(xc, z2, z);
1882                                 }
1883                                 while(x < xae);
1884                         *xc = carry;
1885                         }
1886                 if (y = *xb >> 16) {
1887                         x = xa;
1888                         xc = xc0;
1889                         carry = 0;
1890                         z2 = *xc;
1891                         do {
1892                                 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
1893                                 carry = z >> 16;
1894                                 Storeinc(xc, z, z2);
1895                                 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
1896                                 carry = z2 >> 16;
1897                                 }
1898                                 while(x < xae);
1899                         *xc = z2;
1900                         }
1901                 }
1902 #else
1903         for(; xb < xbe; xc0++) {
1904                 if (y = *xb++) {
1905                         x = xa;
1906                         xc = xc0;
1907                         carry = 0;
1908                         do {
1909                                 z = *x++ * y + *xc + carry;
1910                                 carry = z >> 16;
1911                                 *xc++ = z & 0xffff;
1912                                 }
1913                                 while(x < xae);
1914                         *xc = carry;
1915                         }
1916                 }
1917 #endif
1918 #endif
1919         for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
1920         c->wds = wc;
1921         return c;
1922         }
1923
1924  static Bigint *
1925 pow5mult(Bigint *b, int k MTd)
1926 {
1927         Bigint *b1, *p5, *p51;
1928 #ifdef MULTIPLE_THREADS
1929         ThInfo *TI;
1930 #endif
1931         int i;
1932         static int p05[3] = { 5, 25, 125 };
1933
1934         if ((i = k & 3))
1935                 b = multadd(b, p05[i-1], 0 MTa);
1936
1937         if (!(k >>= 2))
1938                 return b;
1939 #ifdef  MULTIPLE_THREADS
1940         if (!(TI = *PTI))
1941                 *PTI = TI = get_TI();
1942 #endif
1943         if (!(p5 = p5s)) {
1944                 /* first time */
1945 #ifdef MULTIPLE_THREADS
1946                 if (!(TI = *PTI))
1947                         *PTI = TI = get_TI();
1948                 if (TI == &TI0)
1949                         ACQUIRE_DTOA_LOCK(1);
1950                 if (!(p5 = p5s)) {
1951                         p5 = p5s = i2b(625 MTa);
1952                         p5->next = 0;
1953                         }
1954                 if (TI == &TI0)
1955                         FREE_DTOA_LOCK(1);
1956 #else
1957                 p5 = p5s = i2b(625 MTa);
1958                 p5->next = 0;
1959 #endif
1960                 }
1961         for(;;) {
1962                 if (k & 1) {
1963                         b1 = mult(b, p5 MTa);
1964                         Bfree(b MTa);
1965                         b = b1;
1966                         }
1967                 if (!(k >>= 1))
1968                         break;
1969                 if (!(p51 = p5->next)) {
1970 #ifdef MULTIPLE_THREADS
1971                         if (!TI && !(TI = *PTI))
1972                                 *PTI = TI = get_TI();
1973                         if (TI == &TI0)
1974                                 ACQUIRE_DTOA_LOCK(1);
1975                         if (!(p51 = p5->next)) {
1976                                 p51 = p5->next = mult(p5,p5 MTa);
1977                                 p51->next = 0;
1978                                 }
1979                         if (TI == &TI0)
1980                                 FREE_DTOA_LOCK(1);
1981 #else
1982                         p51 = p5->next = mult(p5,p5);
1983                         p51->next = 0;
1984 #endif
1985                         }
1986                 p5 = p51;
1987                 }
1988         return b;
1989         }
1990
1991  static Bigint *
1992 lshift(Bigint *b, int k MTd)
1993 {
1994         int i, k1, n, n1;
1995         Bigint *b1;
1996         ULong *x, *x1, *xe, z;
1997
1998 #ifdef Pack_32
1999         n = k >> 5;
2000 #else
2001         n = k >> 4;
2002 #endif
2003         k1 = b->k;
2004         n1 = n + b->wds + 1;
2005         for(i = b->maxwds; n1 > i; i <<= 1)
2006                 k1++;
2007         b1 = Balloc(k1 MTa);
2008         x1 = b1->x;
2009         for(i = 0; i < n; i++)
2010                 *x1++ = 0;
2011         x = b->x;
2012         xe = x + b->wds;
2013 #ifdef Pack_32
2014         if (k &= 0x1f) {
2015                 k1 = 32 - k;
2016                 z = 0;
2017                 do {
2018                         *x1++ = *x << k | z;
2019                         z = *x++ >> k1;
2020                         }
2021                         while(x < xe);
2022                 if ((*x1 = z))
2023                         ++n1;
2024                 }
2025 #else
2026         if (k &= 0xf) {
2027                 k1 = 16 - k;
2028                 z = 0;
2029                 do {
2030                         *x1++ = *x << k  & 0xffff | z;
2031                         z = *x++ >> k1;
2032                         }
2033                         while(x < xe);
2034                 if (*x1 = z)
2035                         ++n1;
2036                 }
2037 #endif
2038         else do
2039                 *x1++ = *x++;
2040                 while(x < xe);
2041         b1->wds = n1 - 1;
2042         Bfree(b MTa);
2043         return b1;
2044         }
2045
2046  static int
2047 cmp(Bigint *a, Bigint *b)
2048 {
2049         ULong *xa, *xa0, *xb, *xb0;
2050         int i, j;
2051
2052         i = a->wds;
2053         j = b->wds;
2054 #ifdef DEBUG
2055         if (i > 1 && !a->x[i-1])
2056                 Bug("cmp called with a->x[a->wds-1] == 0");
2057         if (j > 1 && !b->x[j-1])
2058                 Bug("cmp called with b->x[b->wds-1] == 0");
2059 #endif
2060         if (i -= j)
2061                 return i;
2062         xa0 = a->x;
2063         xa = xa0 + j;
2064         xb0 = b->x;
2065         xb = xb0 + j;
2066         for(;;) {
2067                 if (*--xa != *--xb)
2068                         return *xa < *xb ? -1 : 1;
2069                 if (xa <= xa0)
2070                         break;
2071                 }
2072         return 0;
2073         }
2074
2075  static Bigint *
2076 diff(Bigint *a, Bigint *b MTd)
2077 {
2078         Bigint *c;
2079         int i, wa, wb;
2080         ULong *xa, *xae, *xb, *xbe, *xc;
2081 #ifdef ULLong
2082         ULLong borrow, y;
2083 #else
2084         ULong borrow, y;
2085 #ifdef Pack_32
2086         ULong z;
2087 #endif
2088 #endif
2089
2090         i = cmp(a,b);
2091         if (!i) {
2092                 c = Balloc(0 MTa);
2093                 c->wds = 1;
2094                 c->x[0] = 0;
2095                 return c;
2096                 }
2097         if (i < 0) {
2098                 c = a;
2099                 a = b;
2100                 b = c;
2101                 i = 1;
2102                 }
2103         else
2104                 i = 0;
2105         c = Balloc(a->k MTa);
2106         c->sign = i;
2107         wa = a->wds;
2108         xa = a->x;
2109         xae = xa + wa;
2110         wb = b->wds;
2111         xb = b->x;
2112         xbe = xb + wb;
2113         xc = c->x;
2114         borrow = 0;
2115 #ifdef ULLong
2116         do {
2117                 y = (ULLong)*xa++ - *xb++ - borrow;
2118                 borrow = y >> 32 & (ULong)1;
2119                 *xc++ = y & FFFFFFFF;
2120                 }
2121                 while(xb < xbe);
2122         while(xa < xae) {
2123                 y = *xa++ - borrow;
2124                 borrow = y >> 32 & (ULong)1;
2125                 *xc++ = y & FFFFFFFF;
2126                 }
2127 #else
2128 #ifdef Pack_32
2129         do {
2130                 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
2131                 borrow = (y & 0x10000) >> 16;
2132                 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
2133                 borrow = (z & 0x10000) >> 16;
2134                 Storeinc(xc, z, y);
2135                 }
2136                 while(xb < xbe);
2137         while(xa < xae) {
2138                 y = (*xa & 0xffff) - borrow;
2139                 borrow = (y & 0x10000) >> 16;
2140                 z = (*xa++ >> 16) - borrow;
2141                 borrow = (z & 0x10000) >> 16;
2142                 Storeinc(xc, z, y);
2143                 }
2144 #else
2145         do {
2146                 y = *xa++ - *xb++ - borrow;
2147                 borrow = (y & 0x10000) >> 16;
2148                 *xc++ = y & 0xffff;
2149                 }
2150                 while(xb < xbe);
2151         while(xa < xae) {
2152                 y = *xa++ - borrow;
2153                 borrow = (y & 0x10000) >> 16;
2154                 *xc++ = y & 0xffff;
2155                 }
2156 #endif
2157 #endif
2158         while(!*--xc)
2159                 wa--;
2160         c->wds = wa;
2161         return c;
2162         }
2163
2164  static double
2165 ulp(U *x)
2166 {
2167         Long L;
2168         U u;
2169
2170         L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
2171 #ifndef Avoid_Underflow
2172 #ifndef Sudden_Underflow
2173         if (L > 0) {
2174 #endif
2175 #endif
2176 #ifdef IBM
2177                 L |= Exp_msk1 >> 4;
2178 #endif
2179                 word0(&u) = L;
2180                 word1(&u) = 0;
2181 #ifndef Avoid_Underflow
2182 #ifndef Sudden_Underflow
2183                 }
2184         else {
2185                 L = -L >> Exp_shift;
2186                 if (L < Exp_shift) {
2187                         word0(&u) = 0x80000 >> L;
2188                         word1(&u) = 0;
2189                         }
2190                 else {
2191                         word0(&u) = 0;
2192                         L -= Exp_shift;
2193                         word1(&u) = L >= 31 ? 1 : 1 << 31 - L;
2194                         }
2195                 }
2196 #endif
2197 #endif
2198         return dval(&u);
2199         }
2200
2201  static double
2202 b2d(Bigint *a, int *e)
2203 {
2204         ULong *xa, *xa0, w, y, z;
2205         int k;
2206         U d;
2207 #ifdef VAX
2208         ULong d0, d1;
2209 #else
2210 #define d0 word0(&d)
2211 #define d1 word1(&d)
2212 #endif
2213
2214         xa0 = a->x;
2215         xa = xa0 + a->wds;
2216         y = *--xa;
2217 #ifdef DEBUG
2218         if (!y) Bug("zero y in b2d");
2219 #endif
2220         k = hi0bits(y);
2221         *e = 32 - k;
2222 #ifdef Pack_32
2223         if (k < Ebits) {
2224                 d0 = Exp_1 | y >> (Ebits - k);
2225                 w = xa > xa0 ? *--xa : 0;
2226                 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
2227                 goto ret_d;
2228                 }
2229         z = xa > xa0 ? *--xa : 0;
2230         if (k -= Ebits) {
2231                 d0 = Exp_1 | y << k | z >> (32 - k);
2232                 y = xa > xa0 ? *--xa : 0;
2233                 d1 = z << k | y >> (32 - k);
2234                 }
2235         else {
2236                 d0 = Exp_1 | y;
2237                 d1 = z;
2238                 }
2239 #else
2240         if (k < Ebits + 16) {
2241                 z = xa > xa0 ? *--xa : 0;
2242                 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
2243                 w = xa > xa0 ? *--xa : 0;
2244                 y = xa > xa0 ? *--xa : 0;
2245                 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
2246                 goto ret_d;
2247                 }
2248         z = xa > xa0 ? *--xa : 0;
2249         w = xa > xa0 ? *--xa : 0;
2250         k -= Ebits + 16;
2251         d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
2252         y = xa > xa0 ? *--xa : 0;
2253         d1 = w << k + 16 | y << k;
2254 #endif
2255  ret_d:
2256 #ifdef VAX
2257         word0(&d) = d0 >> 16 | d0 << 16;
2258         word1(&d) = d1 >> 16 | d1 << 16;
2259 #else
2260 #undef d0
2261 #undef d1
2262 #endif
2263         return dval(&d);
2264         }
2265
2266  static Bigint *
2267 d2b(U *d, int *e, int *bits MTd)
2268 {
2269         Bigint *b;
2270         int de, k;
2271         ULong *x, y, z;
2272 #ifndef Sudden_Underflow
2273         int i;
2274 #endif
2275 #ifdef VAX
2276         ULong d0, d1;
2277         d0 = word0(d) >> 16 | word0(d) << 16;
2278         d1 = word1(d) >> 16 | word1(d) << 16;
2279 #else
2280 #define d0 word0(d)
2281 #define d1 word1(d)
2282 #endif
2283
2284 #ifdef Pack_32
2285         b = Balloc(1 MTa);
2286 #else
2287         b = Balloc(2 MTa);
2288 #endif
2289         x = b->x;
2290
2291         z = d0 & Frac_mask;
2292         d0 &= 0x7fffffff;       /* clear sign bit, which we ignore */
2293 #ifdef Sudden_Underflow
2294         de = (int)(d0 >> Exp_shift);
2295 #ifndef IBM
2296         z |= Exp_msk11;
2297 #endif
2298 #else
2299         if ((de = (int)(d0 >> Exp_shift)))
2300                 z |= Exp_msk1;
2301 #endif
2302 #ifdef Pack_32
2303         if ((y = d1)) {
2304                 if ((k = lo0bits(&y))) {
2305                         x[0] = y | z << (32 - k);
2306                         z >>= k;
2307                         }
2308                 else
2309                         x[0] = y;
2310 #ifndef Sudden_Underflow
2311                 i =
2312 #endif
2313                     b->wds = (x[1] = z) ? 2 : 1;
2314                 }
2315         else {
2316                 k = lo0bits(&z);
2317                 x[0] = z;
2318 #ifndef Sudden_Underflow
2319                 i =
2320 #endif
2321                     b->wds = 1;
2322                 k += 32;
2323                 }
2324 #else
2325         if (y = d1) {
2326                 if (k = lo0bits(&y))
2327                         if (k >= 16) {
2328                                 x[0] = y | z << 32 - k & 0xffff;
2329                                 x[1] = z >> k - 16 & 0xffff;
2330                                 x[2] = z >> k;
2331                                 i = 2;
2332                                 }
2333                         else {
2334                                 x[0] = y & 0xffff;
2335                                 x[1] = y >> 16 | z << 16 - k & 0xffff;
2336                                 x[2] = z >> k & 0xffff;
2337                                 x[3] = z >> k+16;
2338                                 i = 3;
2339                                 }
2340                 else {
2341                         x[0] = y & 0xffff;
2342                         x[1] = y >> 16;
2343                         x[2] = z & 0xffff;
2344                         x[3] = z >> 16;
2345                         i = 3;
2346                         }
2347                 }
2348         else {
2349 #ifdef DEBUG
2350                 if (!z)
2351                         Bug("Zero passed to d2b");
2352 #endif
2353                 k = lo0bits(&z);
2354                 if (k >= 16) {
2355                         x[0] = z;
2356                         i = 0;
2357                         }
2358                 else {
2359                         x[0] = z & 0xffff;
2360                         x[1] = z >> 16;
2361                         i = 1;
2362                         }
2363                 k += 32;
2364                 }
2365         while(!x[i])
2366                 --i;
2367         b->wds = i + 1;
2368 #endif
2369 #ifndef Sudden_Underflow
2370         if (de) {
2371 #endif
2372 #ifdef IBM
2373                 *e = (de - Bias - (P-1) << 2) + k;
2374                 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
2375 #else
2376                 *e = de - Bias - (P-1) + k;
2377                 *bits = P - k;
2378 #endif
2379 #ifndef Sudden_Underflow
2380                 }
2381         else {
2382                 *e = de - Bias - (P-1) + 1 + k;
2383 #ifdef Pack_32
2384                 *bits = 32*i - hi0bits(x[i-1]);
2385 #else
2386                 *bits = (i+2)*16 - hi0bits(x[i]);
2387 #endif
2388                 }
2389 #endif
2390         return b;
2391         }
2392 #undef d0
2393 #undef d1
2394
2395  static double
2396 ratio(Bigint *a, Bigint *b)
2397 {
2398         U da, db;
2399         int k, ka, kb;
2400
2401         dval(&da) = b2d(a, &ka);
2402         dval(&db) = b2d(b, &kb);
2403 #ifdef Pack_32
2404         k = ka - kb + 32*(a->wds - b->wds);
2405 #else
2406         k = ka - kb + 16*(a->wds - b->wds);
2407 #endif
2408 #ifdef IBM
2409         if (k > 0) {
2410                 word0(&da) += (k >> 2)*Exp_msk1;
2411                 if (k &= 3)
2412                         dval(&da) *= 1 << k;
2413                 }
2414         else {
2415                 k = -k;
2416                 word0(&db) += (k >> 2)*Exp_msk1;
2417                 if (k &= 3)
2418                         dval(&db) *= 1 << k;
2419                 }
2420 #else
2421         if (k > 0)
2422                 word0(&da) += k*Exp_msk1;
2423         else {
2424                 k = -k;
2425                 word0(&db) += k*Exp_msk1;
2426                 }
2427 #endif
2428         return dval(&da) / dval(&db);
2429         }
2430
2431  static const double
2432 tens[] = {
2433                 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
2434                 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
2435                 1e20, 1e21, 1e22
2436 #ifdef VAX
2437                 , 1e23, 1e24
2438 #endif
2439                 };
2440
2441  static const double
2442 #ifdef IEEE_Arith
2443 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
2444 static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
2445 #ifdef Avoid_Underflow
2446                 9007199254740992.*9007199254740992.e-256
2447                 /* = 2^106 * 1e-256 */
2448 #else
2449                 1e-256
2450 #endif
2451                 };
2452 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
2453 /* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
2454 #define Scale_Bit 0x10
2455 #define n_bigtens 5
2456 #else
2457 #ifdef IBM
2458 bigtens[] = { 1e16, 1e32, 1e64 };
2459 static const double tinytens[] = { 1e-16, 1e-32, 1e-64 };
2460 #define n_bigtens 3
2461 #else
2462 bigtens[] = { 1e16, 1e32 };
2463 static const double tinytens[] = { 1e-16, 1e-32 };
2464 #define n_bigtens 2
2465 #endif
2466 #endif
2467
2468 #undef Need_Hexdig
2469 #ifdef INFNAN_CHECK
2470 #ifndef No_Hex_NaN
2471 #define Need_Hexdig
2472 #endif
2473 #endif
2474
2475 #ifndef Need_Hexdig
2476 #ifndef NO_HEX_FP
2477 #define Need_Hexdig
2478 #endif
2479 #endif
2480
2481 #ifdef Need_Hexdig /*{*/
2482 #if 0
2483 static unsigned char hexdig[256];
2484
2485  static void
2486 htinit(unsigned char *h, unsigned char *s, int inc)
2487 {
2488         int i, j;
2489         for(i = 0; (j = s[i]) !=0; i++)
2490                 h[j] = i + inc;
2491         }
2492
2493  static void
2494 hexdig_init(void)       /* Use of hexdig_init omitted 20121220 to avoid a */
2495                         /* race condition when multiple threads are used. */
2496 {
2497 #define USC (unsigned char *)
2498         htinit(hexdig, USC "0123456789", 0x10);
2499         htinit(hexdig, USC "abcdef", 0x10 + 10);
2500         htinit(hexdig, USC "ABCDEF", 0x10 + 10);
2501         }
2502 #else
2503 static unsigned char hexdig[256] = {
2504         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2505         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2506         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2507         16,17,18,19,20,21,22,23,24,25,0,0,0,0,0,0,
2508         0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0,
2509         0,0,0,0,0,0,0,0,0,0,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,0,0,0,0,0,0,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         };
2521 #endif
2522 #endif /* } Need_Hexdig */
2523
2524 #ifdef INFNAN_CHECK
2525
2526 #ifndef NAN_WORD0
2527 #define NAN_WORD0 0x7ff80000
2528 #endif
2529
2530 #ifndef NAN_WORD1
2531 #define NAN_WORD1 0
2532 #endif
2533
2534  static int
2535 match(const char **sp, const char *t)
2536 {
2537         int c, d;
2538         const char *s = *sp;
2539
2540         while((d = *t++)) {
2541                 if ((c = *++s) >= 'A' && c <= 'Z')
2542                         c += 'a' - 'A';
2543                 if (c != d)
2544                         return 0;
2545                 }
2546         *sp = s + 1;
2547         return 1;
2548         }
2549
2550 #ifndef No_Hex_NaN
2551  static void
2552 hexnan(U *rvp, const char **sp)
2553 {
2554         ULong c, x[2];
2555         const char *s;
2556         int c1, havedig, udx0, xshift;
2557
2558         /**** if (!hexdig['0']) hexdig_init(); ****/
2559         x[0] = x[1] = 0;
2560         havedig = xshift = 0;
2561         udx0 = 1;
2562         s = *sp;
2563         /* allow optional initial 0x or 0X */
2564         while((c = *(const unsigned char*)(s+1)) && c <= ' ')
2565                 ++s;
2566         if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X'))
2567                 s += 2;
2568         while((c = *(const unsigned char*)++s)) {
2569                 if ((c1 = hexdig[c]))
2570                         c  = c1 & 0xf;
2571                 else if (c <= ' ') {
2572                         if (udx0 && havedig) {
2573                                 udx0 = 0;
2574                                 xshift = 1;
2575                                 }
2576                         continue;
2577                         }
2578 #ifdef GDTOA_NON_PEDANTIC_NANCHECK
2579                 else if (/*(*/ c == ')' && havedig) {
2580                         *sp = s + 1;
2581                         break;
2582                         }
2583                 else
2584                         return; /* invalid form: don't change *sp */
2585 #else
2586                 else {
2587                         do {
2588                                 if (/*(*/ c == ')') {
2589                                         *sp = s + 1;
2590                                         break;
2591                                         }
2592                                 } while((c = *++s));
2593                         break;
2594                         }
2595 #endif
2596                 havedig = 1;
2597                 if (xshift) {
2598                         xshift = 0;
2599                         x[0] = x[1];
2600                         x[1] = 0;
2601                         }
2602                 if (udx0)
2603                         x[0] = (x[0] << 4) | (x[1] >> 28);
2604                 x[1] = (x[1] << 4) | c;
2605                 }
2606         if ((x[0] &= 0xfffff) || x[1]) {
2607                 word0(rvp) = Exp_mask | x[0];
2608                 word1(rvp) = x[1];
2609                 }
2610         }
2611 #endif /*No_Hex_NaN*/
2612 #endif /* INFNAN_CHECK */
2613
2614 #ifdef Pack_32
2615 #define ULbits 32
2616 #define kshift 5
2617 #define kmask 31
2618 #else
2619 #define ULbits 16
2620 #define kshift 4
2621 #define kmask 15
2622 #endif
2623
2624 #if !defined(NO_HEX_FP) || defined(Honor_FLT_ROUNDS) /*{*/
2625  static Bigint *
2626 increment(Bigint *b MTd)
2627 {
2628         ULong *x, *xe;
2629         Bigint *b1;
2630
2631         x = b->x;
2632         xe = x + b->wds;
2633         do {
2634                 if (*x < (ULong)0xffffffffL) {
2635                         ++*x;
2636                         return b;
2637                         }
2638                 *x++ = 0;
2639                 } while(x < xe);
2640         {
2641                 if (b->wds >= b->maxwds) {
2642                         b1 = Balloc(b->k+1 MTa);
2643                         Bcopy(b1,b);
2644                         Bfree(b MTa);
2645                         b = b1;
2646                         }
2647                 b->x[b->wds++] = 1;
2648                 }
2649         return b;
2650         }
2651
2652 #endif /*}*/
2653
2654 #ifndef NO_HEX_FP /*{*/
2655
2656  static void
2657 rshift(Bigint *b, int k)
2658 {
2659         ULong *x, *x1, *xe, y;
2660         int n;
2661
2662         x = x1 = b->x;
2663         n = k >> kshift;
2664         if (n < b->wds) {
2665                 xe = x + b->wds;
2666                 x += n;
2667                 if (k &= kmask) {
2668                         n = 32 - k;
2669                         y = *x++ >> k;
2670                         while(x < xe) {
2671                                 *x1++ = (y | (*x << n)) & 0xffffffff;
2672                                 y = *x++ >> k;
2673                                 }
2674                         if ((*x1 = y) !=0)
2675                                 x1++;
2676                         }
2677                 else
2678                         while(x < xe)
2679                                 *x1++ = *x++;
2680                 }
2681         if ((b->wds = x1 - b->x) == 0)
2682                 b->x[0] = 0;
2683         }
2684
2685  static ULong
2686 any_on(Bigint *b, int k)
2687 {
2688         int n, nwds;
2689         ULong *x, *x0, x1, x2;
2690
2691         x = b->x;
2692         nwds = b->wds;
2693         n = k >> kshift;
2694         if (n > nwds)
2695                 n = nwds;
2696         else if (n < nwds && (k &= kmask)) {
2697                 x1 = x2 = x[n];
2698                 x1 >>= k;
2699                 x1 <<= k;
2700                 if (x1 != x2)
2701                         return 1;
2702                 }
2703         x0 = x;
2704         x += n;
2705         while(x > x0)
2706                 if (*--x)
2707                         return 1;
2708         return 0;
2709         }
2710
2711 enum {  /* rounding values: same as FLT_ROUNDS */
2712         Round_zero = 0,
2713         Round_near = 1,
2714         Round_up = 2,
2715         Round_down = 3
2716         };
2717
2718  void
2719 gethex( const char **sp, U *rvp, int rounding, int sign MTd)
2720 {
2721         Bigint *b;
2722         const unsigned char *decpt, *s0, *s, *s1;
2723         Long e, e1;
2724         ULong L, lostbits, *x;
2725         int big, denorm, esign, havedig, k, n, nbits, up, zret;
2726 #ifdef IBM
2727         int j;
2728 #endif
2729         enum {
2730 #ifdef IEEE_Arith /*{{*/
2731                 emax = 0x7fe - Bias - P + 1,
2732                 emin = Emin - P + 1
2733 #else /*}{*/
2734                 emin = Emin - P,
2735 #ifdef VAX
2736                 emax = 0x7ff - Bias - P + 1
2737 #endif
2738 #ifdef IBM
2739                 emax = 0x7f - Bias - P
2740 #endif
2741 #endif /*}}*/
2742                 };
2743 #ifdef USE_LOCALE
2744         int i;
2745 #ifdef NO_LOCALE_CACHE
2746         const unsigned char *decimalpoint = (unsigned char*)
2747                 localeconv()->decimal_point;
2748 #else
2749         const unsigned char *decimalpoint;
2750         static unsigned char *decimalpoint_cache;
2751         if (!(s0 = decimalpoint_cache)) {
2752                 s0 = (unsigned char*)localeconv()->decimal_point;
2753                 if ((decimalpoint_cache = (unsigned char*)
2754                                 MALLOC(strlen((const char*)s0) + 1))) {
2755                         strcpy((char*)decimalpoint_cache, (const char*)s0);
2756                         s0 = decimalpoint_cache;
2757                         }
2758                 }
2759         decimalpoint = s0;
2760 #endif
2761 #endif
2762
2763         /**** if (!hexdig['0']) hexdig_init(); ****/
2764         havedig = 0;
2765         s0 = *(const unsigned char **)sp + 2;
2766         while(s0[havedig] == '0')
2767                 havedig++;
2768         s0 += havedig;
2769         s = s0;
2770         decpt = 0;
2771         zret = 0;
2772         e = 0;
2773         if (hexdig[*s])
2774                 havedig++;
2775         else {
2776                 zret = 1;
2777 #ifdef USE_LOCALE
2778                 for(i = 0; decimalpoint[i]; ++i) {
2779                         if (s[i] != decimalpoint[i])
2780                                 goto pcheck;
2781                         }
2782                 decpt = s += i;
2783 #else
2784                 if (*s != '.')
2785                         goto pcheck;
2786                 decpt = ++s;
2787 #endif
2788                 if (!hexdig[*s])
2789                         goto pcheck;
2790                 while(*s == '0')
2791                         s++;
2792                 if (hexdig[*s])
2793                         zret = 0;
2794                 havedig = 1;
2795                 s0 = s;
2796                 }
2797         while(hexdig[*s])
2798                 s++;
2799 #ifdef USE_LOCALE
2800         if (*s == *decimalpoint && !decpt) {
2801                 for(i = 1; decimalpoint[i]; ++i) {
2802                         if (s[i] != decimalpoint[i])
2803                                 goto pcheck;
2804                         }
2805                 decpt = s += i;
2806 #else
2807         if (*s == '.' && !decpt) {
2808                 decpt = ++s;
2809 #endif
2810                 while(hexdig[*s])
2811                         s++;
2812                 }/*}*/
2813         if (decpt)
2814                 e = -(((Long)(s-decpt)) << 2);
2815  pcheck:
2816         s1 = s;
2817         big = esign = 0;
2818         switch(*s) {
2819           case 'p':
2820           case 'P':
2821                 switch(*++s) {
2822                   case '-':
2823                         esign = 1;
2824                         /* no break */
2825                   case '+':
2826                         s++;
2827                   }
2828                 if ((n = hexdig[*s]) == 0 || n > 0x19) {
2829                         s = s1;
2830                         break;
2831                         }
2832                 e1 = n - 0x10;
2833                 while((n = hexdig[*++s]) !=0 && n <= 0x19) {
2834                         if (e1 & 0xf8000000)
2835                                 big = 1;
2836                         e1 = 10*e1 + n - 0x10;
2837                         }
2838                 if (esign)
2839                         e1 = -e1;
2840                 e += e1;
2841           }
2842         *sp = (char*)s;
2843         if (!havedig)
2844                 *sp = (char*)s0 - 1;
2845         if (zret)
2846                 goto retz1;
2847         if (big) {
2848                 if (esign) {
2849 #ifdef IEEE_Arith
2850                         switch(rounding) {
2851                           case Round_up:
2852                                 if (sign)
2853                                         break;
2854                                 goto ret_tiny;
2855                           case Round_down:
2856                                 if (!sign)
2857                                         break;
2858                                 goto ret_tiny;
2859                           }
2860 #endif
2861                         goto retz;
2862 #ifdef IEEE_Arith
2863  ret_tinyf:
2864                         Bfree(b MTa);
2865  ret_tiny:
2866                         Set_errno(ERANGE);
2867                         word0(rvp) = 0;
2868                         word1(rvp) = 1;
2869                         return;
2870 #endif /* IEEE_Arith */
2871                         }
2872                 switch(rounding) {
2873                   case Round_near:
2874                         goto ovfl1;
2875                   case Round_up:
2876                         if (!sign)
2877                                 goto ovfl1;
2878                         goto ret_big;
2879                   case Round_down:
2880                         if (sign)
2881                                 goto ovfl1;
2882                         goto ret_big;
2883                   }
2884  ret_big:
2885                 word0(rvp) = Big0;
2886                 word1(rvp) = Big1;
2887                 return;
2888                 }
2889         n = s1 - s0 - 1;
2890         for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
2891                 k++;
2892         b = Balloc(k MTa);
2893         x = b->x;
2894         n = 0;
2895         L = 0;
2896 #ifdef USE_LOCALE
2897         for(i = 0; decimalpoint[i+1]; ++i);
2898 #endif
2899         while(s1 > s0) {
2900 #ifdef USE_LOCALE
2901                 if (*--s1 == decimalpoint[i]) {
2902                         s1 -= i;
2903                         continue;
2904                         }
2905 #else
2906                 if (*--s1 == '.')
2907                         continue;
2908 #endif
2909                 if (n == ULbits) {
2910                         *x++ = L;
2911                         L = 0;
2912                         n = 0;
2913                         }
2914                 L |= (hexdig[*s1] & 0x0f) << n;
2915                 n += 4;
2916                 }
2917         *x++ = L;
2918         b->wds = n = x - b->x;
2919         n = ULbits*n - hi0bits(L);
2920         nbits = Nbits;
2921         lostbits = 0;
2922         x = b->x;
2923         if (n > nbits) {
2924                 n -= nbits;
2925                 if (any_on(b,n)) {
2926                         lostbits = 1;
2927                         k = n - 1;
2928                         if (x[k>>kshift] & 1 << (k & kmask)) {
2929                                 lostbits = 2;
2930                                 if (k > 0 && any_on(b,k))
2931                                         lostbits = 3;
2932                                 }
2933                         }
2934                 rshift(b, n);
2935                 e += n;
2936                 }
2937         else if (n < nbits) {
2938                 n = nbits - n;
2939                 b = lshift(b, n MTa);
2940                 e -= n;
2941                 x = b->x;
2942                 }
2943         if (e > emax) {
2944  ovfl:
2945                 Bfree(b MTa);
2946  ovfl1:
2947                 Set_errno(ERANGE);
2948 #ifdef Honor_FLT_ROUNDS
2949                 switch (rounding) {
2950                   case Round_zero:
2951                         goto ret_big;
2952                   case Round_down:
2953                         if (!sign)
2954                                 goto ret_big;
2955                         break;
2956                   case Round_up:
2957                         if (sign)
2958                                 goto ret_big;
2959                   }
2960 #endif
2961                 word0(rvp) = Exp_mask;
2962                 word1(rvp) = 0;
2963                 return;
2964                 }
2965         denorm = 0;
2966         if (e < emin) {
2967                 denorm = 1;
2968                 n = emin - e;
2969                 if (n >= nbits) {
2970 #ifdef IEEE_Arith /*{*/
2971                         switch (rounding) {
2972                           case Round_near:
2973                                 if (n == nbits && (n < 2 || lostbits || any_on(b,n-1)))
2974                                         goto ret_tinyf;
2975                                 break;
2976                           case Round_up:
2977                                 if (!sign)
2978                                         goto ret_tinyf;
2979                                 break;
2980                           case Round_down:
2981                                 if (sign)
2982                                         goto ret_tinyf;
2983                           }
2984 #endif /* } IEEE_Arith */
2985                         Bfree(b MTa);
2986  retz:
2987                         Set_errno(ERANGE);
2988  retz1:
2989                         rvp->d = 0.;
2990                         return;
2991                         }
2992                 k = n - 1;
2993                 if (lostbits)
2994                         lostbits = 1;
2995                 else if (k > 0)
2996                         lostbits = any_on(b,k);
2997                 if (x[k>>kshift] & 1 << (k & kmask))
2998                         lostbits |= 2;
2999                 nbits -= n;
3000                 rshift(b,n);
3001                 e = emin;
3002                 }
3003         if (lostbits) {
3004                 up = 0;
3005                 switch(rounding) {
3006                   case Round_zero:
3007                         break;
3008                   case Round_near:
3009                         if (lostbits & 2
3010                          && (lostbits & 1) | (x[0] & 1))
3011                                 up = 1;
3012                         break;
3013                   case Round_up:
3014                         up = 1 - sign;
3015                         break;
3016                   case Round_down:
3017                         up = sign;
3018                   }
3019                 if (up) {
3020                         k = b->wds;
3021                         b = increment(b MTa);
3022                         x = b->x;
3023                         if (denorm) {
3024 #if 0
3025                                 if (nbits == Nbits - 1
3026                                  && x[nbits >> kshift] & 1 << (nbits & kmask))
3027                                         denorm = 0; /* not currently used */
3028 #endif
3029                                 }
3030                         else if (b->wds > k
3031                          || ((n = nbits & kmask) !=0
3032                              && hi0bits(x[k-1]) < 32-n)) {
3033                                 rshift(b,1);
3034                                 if (++e > Emax)
3035                                         goto ovfl;
3036                                 }
3037                         }
3038                 }
3039 #ifdef IEEE_Arith
3040         if (denorm)
3041                 word0(rvp) = b->wds > 1 ? b->x[1] & ~0x100000 : 0;
3042         else
3043                 word0(rvp) = (b->x[1] & ~0x100000) | ((e + 0x3ff + 52) << 20);
3044         word1(rvp) = b->x[0];
3045 #endif
3046 #ifdef IBM
3047         if ((j = e & 3)) {
3048                 k = b->x[0] & ((1 << j) - 1);
3049                 rshift(b,j);
3050                 if (k) {
3051                         switch(rounding) {
3052                           case Round_up:
3053                                 if (!sign)
3054                                         increment(b);
3055                                 break;
3056                           case Round_down:
3057                                 if (sign)
3058                                         increment(b);
3059                                 break;
3060                           case Round_near:
3061                                 j = 1 << (j-1);
3062                                 if (k & j && ((k & (j-1)) | lostbits))
3063                                         increment(b);
3064                           }
3065                         }
3066                 }
3067         e >>= 2;
3068         word0(rvp) = b->x[1] | ((e + 65 + 13) << 24);
3069         word1(rvp) = b->x[0];
3070 #endif
3071 #ifdef VAX
3072         /* The next two lines ignore swap of low- and high-order 2 bytes. */
3073         /* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */
3074         /* word1(rvp) = b->x[0]; */
3075         word0(rvp) = ((b->x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b->x[1] << 16);
3076         word1(rvp) = (b->x[0] >> 16) | (b->x[0] << 16);
3077 #endif
3078         Bfree(b MTa);
3079         }
3080 #endif /*!NO_HEX_FP}*/
3081
3082  static int
3083 dshift(Bigint *b, int p2)
3084 {
3085         int rv = hi0bits(b->x[b->wds-1]) - 4;
3086         if (p2 > 0)
3087                 rv -= p2;
3088         return rv & kmask;
3089         }
3090
3091  static int
3092 quorem(Bigint *b, Bigint *S)
3093 {
3094         int n;
3095         ULong *bx, *bxe, q, *sx, *sxe;
3096 #ifdef ULLong
3097         ULLong borrow, carry, y, ys;
3098 #else
3099         ULong borrow, carry, y, ys;
3100 #ifdef Pack_32
3101         ULong si, z, zs;
3102 #endif
3103 #endif
3104
3105         n = S->wds;
3106 #ifdef DEBUG
3107         /*debug*/ if (b->wds > n)
3108         /*debug*/       Bug("oversize b in quorem");
3109 #endif
3110         if (b->wds < n)
3111                 return 0;
3112         sx = S->x;
3113         sxe = sx + --n;
3114         bx = b->x;
3115         bxe = bx + n;
3116         q = *bxe / (*sxe + 1);  /* ensure q <= true quotient */
3117 #ifdef DEBUG
3118 #ifdef NO_STRTOD_BIGCOMP
3119         /*debug*/ if (q > 9)
3120 #else
3121         /* An oversized q is possible when quorem is called from bigcomp and */
3122         /* the input is near, e.g., twice the smallest denormalized number. */
3123         /*debug*/ if (q > 15)
3124 #endif
3125         /*debug*/       Bug("oversized quotient in quorem");
3126 #endif
3127         if (q) {
3128                 borrow = 0;
3129                 carry = 0;
3130                 do {
3131 #ifdef ULLong
3132                         ys = *sx++ * (ULLong)q + carry;
3133                         carry = ys >> 32;
3134                         y = *bx - (ys & FFFFFFFF) - borrow;
3135                         borrow = y >> 32 & (ULong)1;
3136                         *bx++ = y & FFFFFFFF;
3137 #else
3138 #ifdef Pack_32
3139                         si = *sx++;
3140                         ys = (si & 0xffff) * q + carry;
3141                         zs = (si >> 16) * q + (ys >> 16);