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