src/base/third_party/dmg_fp/dtoa.cc

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