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