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