wsutil/dtoa.c

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